The High Speed Flow of Gas Around 
Blunt Bodies 


HYMAN SERBIN 
(Hughes Aircraft Company, Culver City, California) 


SuMMaRY: A number of results derived by the author in earlier reports* on the 

flow of air around blunt bodies moving at high speed are here collected in a 

unified analysis. The theory predicts in a satisfactory way the shock shape 

and detachment distance for two blunt bodies, a flat disc and a sphere. It is 

shown that the density ratio across a normak shock is a useful parameter. 

combining the effects of both the free stream Mach number and the ratio of 
specific heats. 


Introduction 


When a body moves through the air, there is a transfer of momentum and 
nergy from the air to the body. The momentum transfer appears as a drag on the 
body and the energy transfer as heat. When the velocity of the body is high, as 
or meteors, the heat input causes surface melting and the material erodes under 


the aerodynamic stresses. 


A sharp-nosed body is particularly vulnerable to the effects of heating. If x 
epresents the distance from the tip of a conical nose to some cross section, then the 
ateral area of the conical tip is of order x* and the volume of order x°. The heat 
iputt is therefore proportional to x°, so that the temperature rise is proportional to 
*/x*=x~', showing that the temperature would increase indefinitely as one 
approached the tip x=0. Therefore a sharp nose would burn off. 


If a body is to be designed to maintain structural integrity at high speeds, then 
t must be blunt. In addition, the speed must be kept low enough so that the body 
Hoes not melt under aerodynamic heating. Therefore, practical considerations of 
gineering design dictate consideration of blunt, high-drag bodies. 


Such shapes are not new to aerodynamic science; they were studied by aero- 
lynamicists, particularly in the early days of the science. One calls to mind 
arman’s theory of vortex streets and Prandtl’s work on separated boundary layers. 


ese earlier reports were prepared during 1956 under the auspices of the RAND Corporation 
of Santa Monica, California. 


This is based on a constant rate of heat transfer. Actually, the rate increases as x decreases, 
again to the disadvantage of the tip. 
Received January 1958. 
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The practical engineer of those days might have called such studies “ academic,” 
because the trend in low-speed design soon led to the slender, low-drag bodies. 
These bodies in turn proved relatively easy to analyse and the profoundly difficult 
questions related to blunt bodies were brushed aside in the advance of aeronautical 
design. 


It is a curious turn of events to be compelled, now from the design point of 
view, once again to examine the aerodynamics of blunt, high drag bodies, not at 
low speeds but at high speeds. 


The first development in this field goes back to Newton. He proposed a 
theory of aerodynamic resistance in which the gas, regarded as a stream of particles, 
completely independent of one another, makes an inelastic impact on the body. 
The assumptions of the theory are valid when the directed energy of the stream is 
large in comparison with the internal energy; hence the Newtonian theory may be 
expected to be valid in the limiting case of high speed. And indeed, the observed 
pressure distributions on the fore part of blunt bodies do check well with the 
Newtonian values. 


From a gas dynamical point of view, the Newtonian theory may be described 
as one valid for the case y=1, where y is the ratio of specific heats. Also, because 
the body does not transmit a signal ahead into the gas stream, it may be said that 
the shock wave which normally lies ahead of the body is now coincident with the 
fore part of the body. To obtain more information on the flow field behind the 
shock, it is natural to seek a theory in which y is near unity without, however, 
attaining the limit.* 


A second approximation, made simultaneously, is that the free stream Mach 
number M,, is large. It turns out that the theory of the subsonic flow field in the 
limiting cases 

y—> 1, M,,— © 


can be described in terms of a single parameter 


_ density behind normal shock 
~ density in front of normal shock ° 


(1) 


Therefore, to set the theory in proper perspective, all the parameters of the 
subsonic flow field are regarded as dependent on K and it is assumed that, as K 
approaches infinity, the declinations of the streamlines and the shock front (from 
a plane normal to the direction of flow) approach zero like K~'. This assumption 
will be justified @ posteriori by displaying a solution of the complete system of 
equations such that the solution satisfies the equations to a principal order in K 
and has the asymptotic properties already postulated. 


*According to the classical kinetic theory of gases, y=1+2/n, where n is the number of 
degrees of freedom of the molecule. It appears that y is closer to unity the larger n is. a 
condition which is brought about by suitably high temperatures. However, at still higher 
temperatures, complete dissociation takes place, leading to the value 5/3 for y. 
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mic, The theory given here is applicable for a perfect gas. There is now considerable 
dies, literature on the effect of real gas properties on shock characteristics; however, it 
ficult J does not seem to have been settled whether the gas effects can be suitably described 
tical by an appropriate value of y. Therefore, there is still a question to be settled of 

how to apply the results of the present work to the case where there is extensive 
nt of | dissociation across the shock wave. 


ot at 
NOTATION 


od a a speed of sound 
cles, 


f@)=M (9.9) 
m is K density ratio across a normal shock (> 1) 
y be 


rved 
the n_ distance along an arc normal to the streamlines 


M Mach number (regarded as a function of 9, v) 


p Static pressure 

q_ velocity 

R_ radius of disc or sphere, in Sections 5 and 6; gas constant, elsewhere 
r radial distance of a point to the axis of symmetry 

r, radial distance to sonic point on the shock produced by a sphere 
S entropy 

s arc length along a streamline 

T temperature 

x abscissa 

a shock wave angle 

y ratio of specific heats 

4s shock stand-off distance 


6 inclination of the streamline 
6’ =2/2-6 


€=Kix/R 
p density of the gas 
o=K? (r/R) 
® first co-ordinate of a point in the flow field 
y¥ stream function, second co-ordinate of a point in the flow field 
Suffixes 
co denotes conditions ahead of shock 


n,$,%,% denote partial differentiation with respect to the indicated suffix. 


November 1958 315 


ibed 

use 

that 

the 

the 

ver, 

ach 

the 

(1) 

K 

om 

ion 

of 

K 

of 

.a 

her 

erly 


HYMAN SERBIN 


2. Normal Shocks at Hypersonic Speeds 


If the density ratio across a normal shock is specified by the parameter 
K (> 1), then the velocity ratio is given by 


(2) 
The Rankine formula for pressure 
Poo (1 + p (1 + 
becomes, in the limit y —> 1, M,, —> 00, M—>0, 


Using equations (1) and (3), the square of the speed of sound takes the form 
(y—> I) 


(4) 


The entropy rise, AS, across the shock is given in general in terms of the free 
stream Mach number by the formula 


(y—1) +2 


jog [ yR (y+ 1) M,.” ]. 


y-1 y+1 (9) 


In the limit y —> 1, this simplifies to 


2M. 


A change 6M, in free stream Mach number M,, therefore leads to a change 6S in 
entropy behind the shock given by 


RM.,,6M,, 
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— Shock front 


q 


Ficure 1. Oblique shock. 


However, the density ratio K across the shock takes the form, when y —> 1, 


Hence, equation (7) can be expressed in the form 


3. Plane Oblique Shocks at Hypersonic Speeds 


Plane oblique shocks with normals inclined only slightly to the direction of 
flow are now considered. The components of velocity in front and behind the 
shock are shown in Fig. 1 for the approximation <’ small. 


Referring to the triangle behind the shock, there exist the geometric relations 
(10) 


Changing from 2 and 6 to the complements 2’ and @ respectively, and using the 
approximation that 2’ and 6 are small, equation (10) can be put in the form 


1 


(12) 


Now the tentative assumption* is made that 6’ and 2’ vary like K~!. Then, in 
equation (11), the first term on the right hand side is negligible, and hence 


*The results will turn out to be consistent with this assumption. 
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The Mach number behind the shock can now be calculated, using equations (4) 
and (13), 


Solving equations (12) and (14) for 2’ and @, 


Equations (15) and (16) constitute a complete description of the shock boundary 
conditions of the flow field. 


4. The Flow Field in Front of a Circular Disc 


The flow field produccd by a blunt body in a high speed flow, such as that 
shown in Fig. 2, is now considered. Part of the shock front will be inclined almost 
normal to the flow and the density ratio across the part of the shock will be almost 
constant, namely K. The discussion henceforth will be confined to the flow field 
directly behind the steeply inclined shock front. 


It is convenient to write the equations of fluid motion in terms of “ natural 
co-ordinates.” 


At each point 7 of the flow field (Fig. 2), place a rectangular co-ordinate 
system (n, s) such that the origin of the co-ordinate system coincides with point P. 
the positive s-axis coincides with the local velocity vector, and the positive n-axis 
is directed towards the concave side of streamline. Each point near T has co- 
ordinates (m, s) and the velocity q and flow inclination 6 may be regarded as 
functions of m and s. Suffixes are used to denote partial differentiation. 


In these co-ordinates, the equations of motion take the form 


l T 
q 


(M* — 1) cos 6”, . (18) 
q 


The entropy S is constant along a streamline, so that the change in S required to 
form S, may be determined at point S (Fig. 2) where the streamline intersects the 
shock. Equation (9) shows the change in entropy across a normal shock due to a 
change in Mach number. 
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FiGuRE 2. Schematic diagram of physical flow field (cross section). 


For an oblique shock, M., in equation (9) must be replaced by M,, sin z, so 
that equation (9) takes the form 


éS=RK cos a da 


(19) 


All quantities in equation (20), with the exception of 2’, are evaluated at the point T 
in the flow field; z’ is evaluated at the point V (Fig. 2). 


The co-ordinates s and n are local co-ordinates, useful only for points near T. 
A co-ordinate system useful for the whole flow field behind the shock can be intro- 
duced with the help of the stream function ¥ of the flow. ¥y is defined by writing 


v=0 on ORQ (Fig. 2). & 


When dn=0, then dy =0, i.e. ¥ is constant along a streamline such as VT. 
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Now construct the family of curves, such as WT, orthogonal to the streamlines. 
Each such curve will intersect the shock at a point, such as W. Denote the value 
of the stream function at the point V (or T) by ¥. Then the pair (9, ¥) serves to 
iocate the point T. 


The point P on the shock represents the intersection of the sonic line with the 
shock. The unit of stream function will be chosen so that Y=1 at P. 


With the choice of unit for » and y, the three-dimensional physical region 
behind the shock is mapped into a region in the (9, ¥)-plane, as shown in Fig. 3. 
Points in the physical space are represented by the corresponding primed (’) letters 
in the new representation. The wall appears as a part of the ¢-axis and the shock 
as part of the line ¢=¥; the sonic line is represented schematically by the broken 
curve. The transformation is degenerate for points on the axis of symmetry of the 
flow, all such points coinciding with the origin ¢=y~=0. 


o’R’ 


0) 


FiGuRE 3. Schematic diagram of flow region in (9, ¥)-plane. 


All the field variables such as q, 6’, etc., will now be represented by functions 
of (9, ¥). The transformation from (s, n) to (¢, ¥) is simplified when K is large. 
For, in that case, according to the tentative hypothesis introduced in Section 3, 
6’ is small and the streamlines in Fig. 2 are almost vertical. It follows in first 
approximation that 


Furthermore, if the unit of length is chosen as the radial distance of point P from 
the axis of symmetry, then, approximately, 


In moving along a streamline, the gas undergoes an isentropic process. 
Therefore 
q MT}, 
where T oc 
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Ficure 6(a). Comparison of experimental and theoretical shock shapes. 
Circular disc. 
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FIGURE 6(b). Comparison of experimental and theoretical shock shapes. Sphere. (M._=5°8). 
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Ledley _ MY) [ 1+3G-DM*q, 


M ¥) exp { [M? ¥)— M? (, ¥)] } (25) 


Mv, ¥) 


as y—> 


In a similar way, it follows that the temperature is constant along the stream- 
line, and hence, by equation (4), is constant everywhere behind the shock. 
Therefore, behind the shock, 


dq dM 


as y—>1. Across the shock, equations (1), (13), and (15) may be applied to give 


. . . . QD 


Using equations (13), (24), (25), (26), and (27), equation (21) can be put in the form 


dy = (9, OOM dn 


But, in the units chosen, 
TP oo = 1. 


The entropy term, equation (20), can be transformed, using equation (15), to 


g” dy dn’ 


Use of this expression in equation (17), changing from (n, s) to (¢, ¥) with the aid 
of equations (23), (24), and (28), and changing the dependent variable from g to M 
by equation (26), results in the following relation : — 


1 OM (9, ve ow” __ MU, dM(y,v) 


0g M?(9,¥) 
(29) 


Since 6 oc K~-, the second term on the left hand side of equation (29) is of 
order K~* and is therefore negligible in comparison with the other terms in that 
equation. 
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dM (¥, ¥) 


1, OM (9.¥) _ 
Hence M (9, au =M (v, 


which may be integrated from Y=¥ to Y=9, 

M? (», ») — M* (9, ¥)= M? (¢, M? (¥, ¥). 

Hence M (4, =M (y, v), ; (0) 
showing that the Mach number, and therefore also the velocity, are constant along 
each streamline. This might have been expected, in view of the comment already 


made that the temperature is constant everywhere behind the shock, and because 
the flow is isentropic along a streamline. 


The second flow equation, equation (18), then takes the form 


Integrating from Y=0 to Y=4, and noting that @ (9, 0)=0, 


Using equation (16), this result takes the form 


_ 


0 


Denoting M (9, 9) by f (9), transposing the ¢-term, and differentiating, 


This equation has a solution, satisfying the condition f=1 at ¢=1, given by 


poe 


The shape of the shock curve may be obtained by noting that 


324 The Aeronautical Quarterly 


wh: 
Int 
WI 
lin 
If, 
1 tal 
(9, 2 | dy. 
In 
N 
I 
tl 
p 
V 
N 


FLOW AROUND BLUNT BODIES 


Writing x= K equation (33) can be written 


which may be transformed into the differential equation 


af 


Integration on f gives 


t 
1 


0 


When f=1, numerical integration furnishes the value 0-28 for this integral.* 


The theory as developed does not lead to detailed information about the sonic 
line. It is known that 
M (9, 


If, in addition, one supposes that the sonic line occurs near ¢= 1, then equation (28) 
takes the approximate form 


(35) 


Integration from Y=0 to Y=1 gives, at the shoulder, 


1 


0 


Numerical integration furnishes the value 0-75. 


The total stand-off distance A is given by combining the results of equations 
(34) and (35), 
K!:A=0-28 + 0-75 =1-03. ‘ , . (36) 


It should be emphasised that A so obtained is measured in units of the ordinate of 
the sonic point P. The theory does not furnish a relation between this ordinate 
and the radius of the disc. However, with the assumption used here (7 —> 1), it is 
plausible that the sonic point P is near the shoulder, so that the results derived 
would be valid when the radius of the disc is the unit of length. 


*This value coincides with a result obtained by W. Hayes in an informal memorandum 
Sonic Aspects of Hypersonic Flow, Ramo-Wooldridge Corporation, 1955, 
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(K—1)'/2.x/R 
Ficure 4. Theoretical shock wave due to circular disc. 
Since =r’, by equation (24), the shock shape and position are now completely 


described in parametric form (parameter=f=M (y, ¥)) according to equations (32), 
(34), and (36). The shock curve so described is plotted in Fig. 4.* 


5. The Flow Field in Front of a Sphere 


The analysis of the sphere becomes similar to that of the disc when it is 
supposed that, as K —> 00, the radius R of the sphere is allowed to approach 
infinity according to 


All the considerations up to and including equation (31) remain valid. How- 
ever, integration of equation (31) must now take place from the lower limit, 


#,-o=r/R. 


Changing from ¢ to r according to equation (24), equation (31) is then integrated 
to give 


(38) 


*The abscissa shown in Fig. 4 is £ modified by replacing K by K-1 (see Section 6). 
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Differentiating with respect to r, and transforming from r to a new variable c, 
defined by 


(39) 


equation (38) can be written 


df 


which is to be integrated from f=0 to f=1, subject to the initial condition 
o;-,=0. ‘ . (41 


It can be verified directly that a solution of equation (40) and (41) is given 
simply by the equation 


o=xf, ‘ . (42) 


The point (o, f)=(0, 0) is a singular point and it can be verified by the method of 
isoclines that other solutions of equation (40) which pass through this point either 
never attain the Mach number unity or attain it but exhibit a flow reversal. Hence 
the solution, equation (42), is the only one physically acceptable. 


The sonic point on the shock, characterised by f=1, is given at radius r,, 
according to equations (39) and (42), by 


r,j/R=K-. . ‘ . (43) 
Just as for the disc, the shape of the shock is found by integrating 


dé 


dp 7K W=K'f 


=r/R, 


again using equations (39) and (42). This gives 


In particular, when r=r,, this becomes, using equation (43), 
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Again, as for the disc, the distance between the sonic point P (Fig. 2) and the sphere 
is obtained by integrating dn in equation (35) from ¥=0 to Y=1. This gives 


Kin= 5 Je dy 


4 


This result is based on choosing r, as the unit of length; that is, 


Hence = 3K (*) 
which gives, using equation (43), 
n 
R 3K (47) 


The shock detachment distance 4 can now be calculated, using equations (45), 
(47), and (43). 


f=0 


R LR 


The last term on the right hand side is approximated, for small r,/R, 


1 


l 
1— /{1-(r,/R)?} = 9 RP = 


H A_ 


48 

3K" (4) 

Equations (44) and (48) furnish a complete description of the shock shape and 
location. 
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$21.03 (K-I) 2 
(108) 


2 3 4 5 6 
K,DENSITY RATIO ACROSS DETACHED SHOCK 


FiGurE 5(a). Shock detachment distance, Circular disc. 
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FiGure 5(b). Shock detachment distance. Sphere. 
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6. Extension to Supersonic Flow 


The results giving the shock shape and location are expressed in terms of the 
parameter K, which incorporates both the effects of y and M,. The limit con- 
sidered here, y—> 1, M,,—> © corresponds to K—> ©, so that the results 
obtained could be expected to be applicable at most to the case of high density ratio 
across a normal shock. On the other hand, when the shock is weak, the stand-off 
distance should increase, becoming infinite when K=1 (sonic case). The stand-off 
distance’ formulae, equations (36) and (48), show finite distances at this value of K, 
but the situation may be remedied, still within the framework of this analysis, by 
replacing everywhere K by K-—1. Thus the shock detachment distances now 
become 


1:03 (K-—1)-*for the disc . (49) 


= (K-1)-' for the sphere. (50) 


A similar substitution is required for the shock shape. 


The results of equations (49) and (50) are shown plotted against K in Fig. 5. 
A comparison is shown with experimental values obtained by measuring the detach- 
ment distances from schlieren photographs and calculating the density ratio K from 
the test values of M,, and y=1:4. The agreement is generally satisfactory, except 
for one point corresponding to a test on a sphere in a stream of chlorine. [If it is 
supposed that the conditions of this test resulted in excitation of additional degrees 
of molecular freedom, and therefore in a reduction of y, and hence in an increase 
in K, the disagreement would be reduced. 


Another comparison between test and theory is shown in Fig. 6 (pp. 321-2). 
The shock shapes have been calculated for the test Mach number and y=1-4. 
Again the agreement is quite satisfactory. 
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Flexural Vibrations of a Thin- Walled Cylinder 
of Rectangular Cross Section 


E. H. MANSFIELD 


(Royal Aircraft Establishment) 


SUMMARY: The natural transverse vibrations of a long cylindrical box of doubly 

symmetrical rectangular cross section are considered. Explicit stress-function 

solutions are obtained for the webs and the top and bottom surfaces so that 

the effects of shear lag and shear deflection are inherently included. The 

results are expressed simply in terms of an effective flexural rigidity, which may 
be determined with the aid of a number of graphs. 


l. Introduction 


In considering the flexural vibrations of a beam it is well known that the shear 
tigidity of the webs must be taken into account, particularly when estimating the 
higher modes“: ». For a box beam, the effect of shear lag in the top and bottom 
surfaces must also be considered. This problem has been discussed by Anderson 
and Houbolt“, who considered a cantilever box of doubly-symmetrical rectangular 
cross section and simplified the structural characteristics by assuming that the webs 
took only shear and the top and bottom surfaces were represented by sheets, taking 
only shear, with a central equivalent stringer; the cross section was assumed to be 
laintained by closely spaced stiff ribs and, because of this assumption, the precise 
distribution of mass over the cross section was immaterial. This approach has 
recently been refined by Davenport and Kruszewski® and Davenport®’, who have 
assumed that the top and bottom surfaces have multi-stringer characteristics. 


In the present paper only the natural (i.e. sinusoidal) vibrations of a long 
cylindrical box of doubly-symmetrical rectangular cross section are considered and 
the simplification resulting from this restriction enables the structural characteristics 
of the cylinder to be more accurately treated. 


NoTaTION (see Fig. 1) 
A,C constants 
v vertical displacement of booms 
t time 


frequency 
Received May 1958. 
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Nor 


* 
Ficure 1. Notation. 


Structure Properties 


w width of cross section of box 

h_ height of cross section of box 

t thickness of skin of top and bottom surfaces 

ty thickness of skin of vertical surfaces (webs) 

S relative stiffness of stringers, =(stringer area)/(t x stringer pitch) 
R_ relative stiffness of rib booms, =(rib boom area)/(t x rib pitch) 
F section area of longitudinal corner member (spar boom) 
F, contribution of skin and stringers to each boom (see equation (5)) 
E Young’s modulus 
El flexural rigidity of box 

vy Poisson’s ratio 

m mass per unit length of box 

A=} (complete wavelength of vibration) 


Axes 


Ox,Oy Cartesian co-ordinates, Ox measured longitudinally, Oy measured F 4p 
in plane of walls of box from mid-point of wall nm 
= A 
n=ry/X 
Stresses 3 
Tr, Ty, Try Sheet stresses 
Jz,%, direct stress resultants in reinforced sheet 
Stringer stress 
¢ Aijry stress function 
The Aeronautical Quarterly 
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VIBRATIONS OF A CYLINDRICAL BOX 


Non-Dimensional Parameters 
k=1+S5+R+SR (1—- 9") 


a=1+S(1-v’) 
y= 1 +(1 +v){S+R+SR —v*)} 
e=1+R(l1-v’) 


n,= 


n,= 


p= 

p=th/(2A) 

B=2 (F+F,)/(htw) 

I" relative stiffness (see equation (15)) 


Ui. Ves ¥, functions introduced in equations (4) and (7) 


2. Assumptions 
The following assumptions are made : — 


(i) stress-strain relations are linear, 
(ii) buckling does not take place, 


(iii) ribs do not offer any resistance to warping out of their plane and the 
rib webs do not carry direct stresses, 


(iv) the stiffening effect of stringers and rib booms is adequately represented 
by assuming them to be spread out into elastic sheets with equivalent 
directional properties, 

(v) the flexural rigidity of the spar booms is negligible, 

(vi) the mass of the cross section is assumed to be concentrated along the 
four corner flanges (spar booms). 

Assumptions (i) to (v) are standard practice. Assumption (vi) implies that 
there will be no inertia forces applied to the rib webs, so that the rib webs carry 
no load. The influence of rotatory or longitudinal inertia is considered in 
Appendix III. 


3. Effective Stiffness in Flexural Vibration 


3.1. GENERAL 


For an infinitely long beam in which shear lag and shear deflection effects are 
neglected, the frequency of vibration in a sinusoidal mode of half wavelength A is 


given by 
a=(%) {= 
333 
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When shear lag and shear deflection effects are taken into account, it will be 
found convenient to express the results in the form of an effective flexural rigidity 
such that the frequency is to be determined from the equation 


(2) 


The derivation of explicit stress function solutions for the webs and top and 
bottom surfaces is given in Appendix I. From these solutions the frequency 
corresponding to a given wavelength of vibration may be calculated. The results 
are most easily expressed in the form of a reduced or effective flexural rigidity of 
the box. Thus 


hty 
6. 


(El) effective =Eh* F+F, (3) 
1+ +(F = Ys 
) 2 
where = 
3 (sinh p cosh p—p) . 4 


2p* cosh?p 


= (1 sech? p + (1+) p tanh p 
and p =th/(2d). 


The term F, may be regarded as the effective contribution of the skin and 
stringers, and is given by the equation 


(n, tanh n,u—n, tanh n,n) 


(5) 
where the non-dimensional parameters n,, n., 2, y., 2, « and » are defined in the 
Notation. 


For very long wavelengths, the term p tends to zero and the functions y, and ¥, 
tend to unity and y, tends to zero, so that shear deflection effects vanish; the term 
F, tends to wtk/(2=), so that shear lag effects vanish and 


~ (F + S)+ 


which is in agreement with elementary beam theory. 
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os CANN: 


Ficure 2. Variation of ¥, with § and A/w. 


3.2. EFFECT OF SHEAR LAG 
The effect of shear lag may be represented by the ratio 
, effective contribution of skin and stringers 
contribution of skin and stringers under pure bending 
2eF, 


(n,? — N,”) 


n, tanh n,u—n, tanh nop 


from equation (5). 


The function ¥, does not vary greatly with R and it is plotted in Fig. 2 for 
varying values of S and A/w, assuming that R is zero and v=0°3. 


If the skin is not reinforced by stringers, equation (7) reduces to 


tanh sech* 


3.3. EFFECT OF SHEAR DEFLECTION 


The effect of shear deflection in reducing the effective flexural rigidity of the 
box is represented primarily by the factor 


1+ 


which occurs in equation (3). Because of the presence of the term F, in this factor, 
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FiGcure 3. Variation of ¥, and y,’ with A/h. 


shear deflection effects cannot be entirely separated from shear lag effects. Further- 
more, because of the magnitude of the term (F + F,)/(ht.), quite modest values of ¥; 
may imply an appreciable effect due to shear deflection. The function y, is plotted 
in Fig. 3 for varying values of A/h. 
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d/h 


Ficure 4. Variation of y, and y, with A/h. 


The functions ¥, and wv, are plotted in Fig. 4; they become significantly 
different from unity only if \<(6h. The functions arise partly from the vertical 
straining in the spar webs, an effect that has been artificially magnified in impor- 
tance by assuming the mass of the box to be concentrated at the corners. The case 
in which there are assumed to be no vertical strains in the spar webs is considered 
in Section 3.5. 


3.4. SHORT WAVELENGTHS 
For very short wavelengths, the term p tends to infinity and, from equation (3), 
Et, 
(ED) cstective +y) v) ( 


8xEt, 


The frequency is thus independent of h and the mode of vibration corresponds 
to that of a flexible but inextensible line load of intensity 4m along the edge of a 
semi-infinite sheet of thickness 1,. 


3.5. ASSUMPTION OF No VERTICAL STRAIN IN THE SPAR WEBS 


This case is treated because it corresponds more closely to the mass being 
distributed over the spar webs, rather than concentrated along the spar booms. 
The “ closely spaced stiff ribs ” case can also be derived from it. 

It is shown in Appendix IT that if there are no vertical strains in the spar webs 
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(F+F,)¥,'+ 


1+ ‘) 


(ED) ettective =Eh’ 


,__ tanh p’ 


3 , 
(1 v2) (p’)" (p tanh p’) 


wv,’ =2 (1—?) p’ tanh p’ 


and p= 


The functions v,’, v.’ and wv,’ are plotted in Figs. 3 and 5. As would be 
expected from physical considerations, there is little difference between v,’, y.’, v,’, 
and v,, ¥., ¥,, except for the very short wavelengths. What difference there is for 
moderate values of A/w is mainly due to Poisson’s ratio effects. As p tends to 


infinity the equation corresponding to equation (9) in Section 3.4 now becomes 
MhEt,, 
(1+) 
Ehtw 
m(1+y) 


which agrees with the value obtained by selieg theory when the only mode 
of vibration is due to shear deflection. 


(EI) effective —> 


so that (12) 


Ficure 5. Variation of ¥,’ and v,’ with A/h. 
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RELATIVE STIFFNESS 


16 
FiGurE 6. Variation of relative stiffness with \/h for two examples. 


4. Example 


Compare the flexural vibration characteristics of the parallel boxes whose 
cross sections are shown in Fig. 6. The two-web box is specified by 


w=4h, R=¢ 

S=1, 
while in the five-web box the web thicknesses are }f, t, t, t, 4t, so that the total web 
thickness is the same as for the two-web box. Other properties are as for the two- 
web box, so that, considering only a single cell, we have 


R=0, 


For the two-web box, from equations (3) and (7), 


(ED) ettective = Eh? 


so that, from equation (13), the relative stiffness is given by equation (15):— 
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~ 13L 1424, 


(15) 


This has been plotted against the wavelength parameter A/h in Fig. 6. The 
corresponding values of I" for the five-web box are also shown in Fig. 6. It will be 
seen that the five-web box is appreciably stiffer, due to the reduced effects of shear 
lag. The graphs for the values of I on the assumption of no vertical strain in the 
webs, discussed in Section 3.5, are indistinguishable from the graphs in Fig. 6. 


Appendix I 


ExPLiciIT STRESS FUNCTION SOLUTION 


In considering the behaviour of the top or bottom surface a solution will first be 
found when the boundary conditions along the edges are 


‘ (16) 
and o,=sin (tx/X) 


so that 2A is the complete wavelength of vibration, and the maximum value of c, is 
unity. The stress function for the top and bottom surfaces must satisfy the equation 


0*o 


and the stress resultants (defined here as the resultant forces in the stiffened sheet per 
unit length divided by #) are then given by 


and the stringer stress is given by 


The stress function satisfying equations (16) and (17) is 


sing, (20) 


id ( =) (= cosh n,y — cosh n,p cosh 


(n,? —n,?) cosh n,u cosh 
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Here 7, and 7, are the positive roots of the equation 
a—2yn*+en*=0 


The longitudinal stress resultant is given by 


= sin€, . 22 
(n,? — n,?) cosh n,n cosh (22) 


K (= cosh n,u cosh n,» —n.” cosh n,u cosh 


which may be integrated across the width of the surface to give the total direct load 
in the surface, which is 


wtk (n, tanh n,u—n, tanh n,n) sin € (23) 
(n,?- n,,”) 


Now the total load in the top or bottom booms is 2F sin ¢, so that, on comparing with 
expression (23), it will be seen that the skin and stringers contribute to each boom an 
effective increase in area given by 


_ (n, tanh n,u—n, tanh 
2eu (n,* —n,°) 


F, (24) 


If the skin is unreinforced, n,—>n,—>1 and equation (24) becomes 
wt 
F,= — (tanh sech?u) . 
4u : 


Stress Distribution in the Spar Webs 


If the spar webs are assumed to be of plain sheet, the stress function, which varies 
sinusoidally in the x-direction and which is an odd function of y, is given by 


2 
o= (sinh » + An cosh yn) sing, 


where A and C are constants to be determined from the boundary conditions. If the 


spar webs are reinforced, the stress function corresponding to that given by equation 
(26) is 


2 
o= sinh n,’ +A, sinh n,’y)sin€. . ‘ - 


Practically, however, the influence of web reinforcement is small and is neglected 
here. The limiting case of no vertical strain in the webs is treated in Appendix II. 
Along the top edge, the stresses given by equation (26) are 

Oz, e=C{(1+24A) sinh p + Ap cosh p}sin £ | 
Oy, e= —C (sinh p+ Ap cosh p) sin £ 


Try, e= —C{(1 + A) cosh p + Ap sinh p} cos £ 
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and the condition of horizontal equilibrium between web and boom +skin and stringer 
contribution is 


which can be solved to give 


— {cosh p+ Bp (1 +) sinh p} (30) 
~ cosh p +p sinh p + sinh p + Bp? (1+) cosh p 


The constant C is to be eliminated from the vertical equilibrium condition between 
the inertia loads and t,o, .. The vertical displacement of the booms is first needed and 
this will be found from integrating the stress-strain relations 


E — VO'ys 


ox 


dv 
E—=0,- 


dy 


4 2 (149) 


whence 


Ve=- (1+¥—A (1—9)} cosh p+ Ap (I+»)sinh p } sin (31) 


If m is the mass per unit length of the whole box and © is the frequency, v, can be 
written 


Ve=Vmax Sin € sin Nt 


where “tt” here denotes time, so that 


m d*v, 
sin Ot= — 


Vnax Sin € sin Nt 


and the term C sin sin Mt cancels out, leaving, an equation for ?, namely 


Q?— 4nE { ___ sinh p+ Ap cosh p 
{1-+¥ — A(1—»)} cosh p+ Ap (1 +9) sinh p 


where A is given by equation (30). 
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Equation (33) may be expressed more conveniently by stating that the box has 
an effective flexural rigidity given by 


F+F)+y, (Me) 
(ED cttective = Eh?| F+F,\ ta 


we 


p 


2 p® cosh? p 
ia 


cosh?p - +(1 +) p tanh p. 


Appendix II 


WEBS WITH No VERTICAL STRAIN 


This case is treated because it corresponds more closely to a distributed mass along 
the spar webs, and the “closely spaced rigid ribs” case can also be derived from it. 
The analysis can proceed from equation (27) in which the limit is taken as the vertical 
reinforcement tends to infinity, so that 


1/2 
n’—>( ;) 


and n,’—>0, 


and equation (27) becomes 
2 1/2 
o= { A, sinh ( ) n+ ¢ sin €, 


which yields stress resultants 


4 1/2 
721) sin sinh ( ) 
2 
Fy= —sin€ {4, sinh n+ A*n } 


= Tey = ~cosé{ > cosh( n+A, 


Now the horizontal strain is given by 


Edu/dx=(1—v?) 


so that the condition of horizontal equilibrium between web and boom-+skin and 
stringer contribution is 
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whence 
\1/2 > 1/2 
A,+A, -) cosh p+2Bp (1+) sinh p } =0. (39) 


The vertical displacement of the booms may be found from integrating the stress- 
strain relations 


du, dv\_, 
E ay * = 


The equation of vertical equilibrium between the inertia loads and the forces in 
the webs, corresponding to equation (32), is now 


mAQ? A, 
27E 


whence (42) 


and the results given in Section 3.5 follow from equations (2), (39) and (42). 


Closely Spaced Rigid Ribs 


This case can be derived from the preceding results, taking R=oo in the 
expression for F,. 


Now, as R—> 00 


2{1 


a— E and n,—> 0, 


so that wi{l+S *)} HY. 


2(i- 
Appendix Hl 


INFLUENCE OF LONGITUDINAL INERTIA 


The influence of longitudinal inertia may be determined as follows. From Appendix 
I, corresponding to equation (31), 


AC 
- =p +24 sinh p + Ap (1+ cosh p}cos (44) 
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and the condition of horizontal equilibrium at a boom now becomes 


0? 


Equation (45) could have been derived from equation (29) by replacing 


m 
roy 
ie. by { F+ « 


and the true value of J srectiy, is found from equation (3) with F replaced by expression 
(46). The resulting equation may be rearranged and solved as a quadratic, but for 
practical purposes a single iteration is generally adequate. 
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A Rapid Method for Calculating the 
“Off-Design” Performance of Compressors 
and Turbines 


J. H. HORLOCK 


(University Engineering Laboratory, Cambridge)* 


SUMMARY: A model for the calculation of the “off-design” performance of 
turbo-machines has been suggested by Mellor’, who proposed that the machine 
be split into an infinite number of stages, each producing a small change in 
stagnation enthalpy. This note suggests that the overall performance of a 
machine may be obtained in an integral form which permits a rapid calculation 
of performance maps for rotational speeds and flow rates not greatly different 
from design values. The method may be of use in preliminary design studies 
of multi-stage compressors and turbines. 


Analytical Model 


For a compressor or turbine stage in which the outlet angles relative to stator 
and rotor remain constant, independent of the entry flow coefficient, the increase in 
stagnation enthalpy is 


Ah, = Cag, = UAc, 


=U (c,, tan z,—C,, tan z,) 


[U —c, (tan 2, —tan f,)] 


where ¢=c,/U, the flow coefficient, and t=tan z, — tan 6, =constant. 


*Now at the Mechanical Engineering Department, University of Liverpool. 
Received May 1958. 
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“OFF-DESIGN” PERFORMANCE OF TURBO-MACHINES 


If there is little change in kinetic energy from entry to exit, then 


AT, =~ AT and 


At the design point (wa, $4), 


At operating points off-design, 


From (2) and (3), 


va va Vi 


Figure 1 shows typical compressor and turbine stages for which this relation 
is true, whether z,, 8, are positive or negative. For a compressor stage AT, and v, 
are positive; for the turbine stage AT, and y, are negative. 


It is important to note that the off-design performance of such an ideal stage 
is controlled by ya, the stage loading at design. Plots of ¥/v against #/¢,, for 
compressor and turbine stages of varying ¥4, are shown in Fig. 2. 


An equivalent stage is to be considered across which the temperature rise is to 
be infinitesimally small, but the nature of the off-design performance must not be 
changed, i.e. Ya must remain the same. The actual machine may be replaced by a 
machine which has an infinite number of small stages of infinitesimally small axial 
length; the geometry of the blading of one of these small stages is the same as that 
of the real stage at the same axial location (i.e. f and ¥, are the same). To obtain 
an infinitesimally small temperature rise or drop across each stage, the new ideal 
machine must rotate very slowly. 


Thus if the stage in the real machine, producing a temperature rise AT, is 
replaced by m “ infinitesimal ” stages, each producing a temperature rise AT’, 
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TYPICAL COMPRESSOR STAGE 


-ve Bp- ve 


a, +ve a, + ve 


TYPICAL TURBINE STAGE 


+ve ANGLE 
+ve Cy 


STATOR 


FIGURE 1. 


Since ¢ is also the same, 


Se. 

As m—> ©, c,’, U’—> 0; for the same mass flow to be passed, mm compressors’ 
must be used in parallel. 


To summarise, the actual machine must be replaced by a large number of 
compressors (or turbines) in parallel, each rotating very slowly, but each producing 
the same temperature rise as the original machine in an infinite number of stages. 
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FicureE 2. Stage temperature rise characteristics. 


With such a model, equation (4) may be written 


dT 1 ¢ 
(5 dT, ( Wa 


where Ty, Ya, ¢a, Ua are distributions through the machine at the design point, and 
T, », U are the distributions through the machine “ off-design.” 


NOTATION 

Stagnation enthalpy 
specific heat at constant pressure 
stagnation temperature 
static temperature 
density 
pressure 

@,,%, air angles (absolute) 

B,,8, air angles (relative) 
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t=tan 2,—tan 
c, axial velocity 
C. tangential velocity 

U_ blade speed 
Y=C,AT,/U’, stage loading 
¢=c,/U, flow coefficient 

A_ area 

M_ mass flow rate 


N_ engine speed (r.p.m.) 


£=M/M, | 
A=U,/U constants 
A.B,b 


nm» polytropic, or small stage, efficiency 


k ratio of specific heats 


n= 


u=T/Ty=1+¢ 
r pressure ratio 
a) 


I entry to machines 


II exit from machines 


d design point, condition of maximum efficiency. 


This notation is used in Fig. 1, which shows velocity triangles for a compressor 
stage and a turbine stage. 


2. Analysis of Off-Design Performance 


At the design point 


CraA apa =CrA e e (6) 


For a machine operating off-design, but with the same entry conditions 


(py =pr, 


(7) 
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From (6) and (7), 


where M denotes mass flow rate. 


If », is the polytropic efficiency, then 


=(#) 


where n= (* —1. 


n» is the polytropic efficiency, defined in terms of static states, and the positive 
sign is taken for compressors, the negative sign is taken for turbines. If the poly- 
tropic efficiency does not change off-design, then 


From (8) and (11), 


where £=M/Mag. 
From equations (5) and (12), 


From (9) and (10), 


Writing A=U4/U, 


Thus the temperature distribution through the machine (7, the dependent 
variable) is expressed in terms of the temperature distribution through the machine 
at the design point (74, the independent variable). In general ya, the mean stage 
loading, will vary through the machine (¥.=f (T.)). 
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Writing u= T/T , and differentiating, 


dT du 
dT, ° ° (14) 
From equations (13) and (14) 
du B 
h 
where = 
va—1\ é 
—— )= 
A 


3. Solution of the Differential Equation 
A solution of the differential equation is obtained for a machine in which ¥;, 


is constant. For a specified off-design speed and flow rate, the parameters A 
and B are constant and equation (15) becomes 


at. u"du 


In general the ratio 1«=7T/T, will not be greatly different from unity, for off- 


design points not far from design, so u=1+¢, where < is small 


Equation (16) may be integrated to give 
log. = constant 
Id 


ell 


e,=0 
Ney nB+1 nA —(n+1) 7 
~ nA—(n+1) log. [1+ ( A+B-1 ex | 


where suffix “ II” denotes the exit condition. 


Thus if the design temperature ratio is known, equation (17) may be solved 
for the off-design temperature ratio 


Ty Tr 
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The pressure ratio is given by 


n+1 


Another exact solution may be obtained for a machine in which n=2 and 
1—b? where b=constant. These conditions are met approximately 
for an air compressor of stage efficiency »,=6/7, in which the root loading, Wajroot), 
is constant through the machine. The solution to equation (15) is then 


2 E (GA) b log Vv (EA) b + Onin) |= Tra (20) 


Tn (2 (€) b= b—(Tya/ Tra) 


4. Calculations 


Two sets of calculations have been made, one for a compressor and one for 
a turbine. 


Stage-by-stage calculations for the Armstrong Siddeley “ Mamba ” compressor 
showed that ¥, was approximately constant (within +3 per cent) at a flow rate of 
17:8 Ib./sec. and a rotational speed of 14,500 r.p.m. Using a value of ¥s=0-339 
and a polytropic efficiency of 90:6 per cent (n=2-17) in equations (17), (18) and (19), 
the characteristic curves of pressure rise and temperature rise as functions of mass 
flow shown in Fig. 3 were obtained. 


Agreement between the experimental and calculated characteristics is good, 
except for the lowest speed shown. 


Howell and Bonham” have given a correlation of the relative stage efficiency 
Mp/Npa against F=(/¥Y,)/(¢/o.), where suffix “d” now denotes the condition of 
maximum efficiency. Howell’s correlation is given in Fig. 4. If the point of the 
maximum efficiency is taken as the design point considered previously, then, from 
equation (4), 


_ _ (1-%a 
= . (A+B) 


(1 + neq) 


At exit, Pa= 


(23) 


-[A (1 + NEq)+ 
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MASS FLOW 
FiGuRE 3. Mamba compressor characteristics. 
An arithmetic mean value of F through the machine 
Fn — 4 (1 + +B| (24) 
2 
Using the results of the first set of calculations of the “‘ Mamba ” performance 
(assuming constant polytropic efficiency), new mean values of the poly- 
tropic (or small stage) efficiency were found, from equation (24) and Fig. 4. In 
re-calculating the characteristic curves of pressure rise against mass flow (Fig. 5). 
using the new polytropic efficiencies, it was assumed that the temperature rise 
through the machine was the same as in the first set of calculations, but the value 
of n in equation (19) was modified. 
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FiGureE 4. Variation of stage efficiency “off-design” (after Howell). 


It is probable that the remaining differences between theory and experiment 
shown in Fig. 5 are due to inaccurate calculation of the overall temperature rise, 
especially at lower speeds. In practice the variation of ¥/¥4 with ¢/, will not be 
linear because of the increase in outlet angles z,, 8, with decreasing flow coefficient. 


The off-design performance of the compressor may be illustrated further by a 
study of the stage characteristics. The range of operation off-design is easily 
calculated. 


Cr Uy M U4 


4 U Cy M, U 


For the first stage, 


@ Tu \* 
For the rear stage, — =£A ( =£€A (1+ neq). 
Pa Tra 


n 
For the middle stage, (i +5 
Lines of constant ¢/¢4 for the front and rear stages have been plotted in 
Fig. 3, and the operating ranges of the front, middle and rear stages are shown in 


Fig. 6, corresponding to the overall characteristics shown in Fig. 3 for N=12,500, 
14,500 r.p.m. 
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MODIFIED CALCULATION 
EXPERIMENTAL 
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Ficure 5. Mamba compressor characteristics. 
FRONT STAGE MIDDLE STAGE REAR STAGE 


$a 
FiGurE 6. Compressor stage characteristics (corresponding to overall characteristics of Fig. 3). 
Arrows indicate decreasing entry flow. 
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“OFF-DESIGN” PERFORMANCE OF TURBO-MACHINES 
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PRESSURE RATIO 


FicureE 7. Tubine characteristics ¥,= —1, »,=87°5 per cent. 


The front stage operates at lower values of ¢/; and probably approaches 
stall at the lower speed, as the flow is decreased. As the flow rate is increased, the 
middle and rear stages move away from stall, at both speeds. 


The turbine map of Fig. 7 was calculated for a turbine in which Y,= —1 and 
p=87-5 per cent (n=3). Most remarkable is the virtual coincidence of the two 
constant speed curves shown and the flattening of the curve of mass flow against 
pressure ratio above the design point. It should be emphasised that this is not 
“choking ”, for the effect of Mach number is not included in the analysis. Fig. 8 
shows a comparison of design speed characteristics of two turbines, one in which 
Ya=—1, and the other in which ¥,=-—2. Also shown in Figs. 7 and 8 is the 
“Stodola ellipse” relation® between mass flow and pressure ratio, derived in 
the Appendix. 


5. Conclusions 


An analytical method has been given for the calculation of the off-design 
performance of compressors and turbines. The method is rapid (the performance 
curves of Fig. 3 were calculated by hand in one hour), and gives good agreement 
with experiment near the design point. It is less accurate in regions where stalling 
of individual blade rows may occur. 
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FiGuRE 8. Turbine characteristics. 


Appendix 
THE “ ELLIPSE LAW ” FoR A MULTI-STAGE TURBINE 


The “ellipse law” for the relationship between pressure ratio and mass flow of a 
multi-stage turbine may be derived if it is assumed that the relation between ¥//w and 
9/4 is parabolic, 


(25) 
dT (<: “(es 
dT, Col” 


This assumption implies that the multi-stage turbine may be treated as a series 
of adiabatic nozzle flows. 
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“OFF-DESIGN” PERFORMANCE OF TURBO-MACHINES 


From the continuity equation, at design and off-design, 
CraPada=C, rapraA Id» 
C,pA =C,1p;Ay, 


M 
and Pe Set Pa. 


If the polytropic efficiency remains constant off-design, then 


p 


= constant, 


77 (=) = = 3) Gay 


Upon integration, 


M\2 
(7) T,°"+1+constant (C) . (30) 
d 


When 


T,=T, 
M 2 
2n+1 
and M (31) 
ra 
ind Since zt =constant, 


p 


Writing the mass flow in terms of pressure ratios r=p,/p, ra=P,/ Da, 


2n+1 | 
M ain 
(1/r,) "*! | 


S Table I shows that the value of (2n + 1)/(n+1) is little less than 2 for turbines of 
high efficiency, and the relation between M/M, and r is approximately elliptical. 
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TABLE I 


(per cent) n (2n+1)/(n+1) 
1:33 100 3 1-75 
1:33 87:5 3:56 1-775 
1-4 100 2-5 1-72 


87°5 3 


In an isothermal turbine, equations (26), (29) and (33) would become 


Ap pu 
Apa p Coa 
Ap _ ( Pa _ ( Pa (35) 
Ap, M, M, p 

1 

M _ / | (36) 
M, \ | 

ra 
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Matrix Analysis of Shell Structures with 
Flexible Frames 


J. S§. PRZEMIENIECKI, Ph.D., Grad.R.Ae.S. 


(Bristol Aircraft Limited) 


SUMMARY: A simple matrix method is presented for the deflection and stress 
analysis of cylindrical shell structures of arbitrary cross section stiffened by 
flexible frames. The method is an extension to fuselage structures of the 
Matrix Force Method developed by Argyris, in which the internal load system 
in the structure consists of two parts:—(a) synthetic load distribution, 
represented by the matrix b,, satisfying the external and internal equations of 
equilibrium, and (b) self-equilibrating load systems, represented by the matrix 
b,, which are introduced to satisfy compatibility conditions. The magnitudes 
of these self-equilibrating load systems are determined from the generalised 
compatibility equations formulated using the flexibility matrix f for the un- 
assembled elements of the structure. The self-equilibrating systems are non- 
orthogonal, but are arranged in such a way that the mixing between one system 
and another is kept to a minimum and, consequently, the resulting compati- 
bility equations are well-conditioned. The three basic matrices, b,, b,, and f, 
are compiled using only very simple formulae. The matrices b, and b, depend 
on the geometry of the structure, while the flexibility matrix f is a function: of 
geometry and elastic properties. The present analysis is applied to cut-out 
problems in fuselage structures. It can also be used for problems involving 
thermal loading and diffusion of loads in curved panels stiffened by 
flexible frames. 


1. Introduction 


The introduction of high-speed computing machines has revolutionised the 
routine methods hitherto used for stress and stiffness analysis of aircraft structures 
and has made possible large scale numerical calculations for complex and irregular 
structures. For such calculations, however, structural problems must be presented 
in matrix algebra suitable for automatic computation on electronic digital computers 
and this usually calls for special matrix methods adapted to the type of structure 
considered. There is already a large number of papers which have led to the 
present widespread use of matrix methods in structural problems’. In particular 
the work by Langefors®: » on the use of matrix methods in aircraft stress analysis 
should be mentioned. Recently a comprehensive study of two dual matrix methods 
based on either displacements (Displacement Method) or forces (Force Method) 
has been presented by Argyris®. It appears, however, that at least with the present 
types of construction the Force Method is to be preferred for aircraft structures. 
In the Force Method a simple synthetic stress system—referred to as the basic 
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system—in equilibrium with the applied loading is selected, and then a number of 
self-equilibrating stress systems is introduced so that when the basic system is 
combined with the self-equilibrating stress systems the conditions of equilibrium 
and compatibility are satisfied. These self-equilibrating stress systems may, there- 
fore, be regarded as perturbations of a basic system that is in equilibrium with the 
external loading but does not satisfy the compatibility conditions. 


In the present paper an attempt has been made to extend the Matrix Force 
Method, which was developed in detail for wing structures in Ref. 9, to cylindrical 
shell structures of arbitrary cross section stiffened by flexible frames. All external 
loads are equilibrated by any statically admissible stress distribution, represented by 
the matrix b,, and naturally the simplest possible stress distribution is used. Such 
a stress distribution is in equilibrium with the applied loads and it also satisfies 
internal equilibrium conditions, but it violates the compatibility conditions for 
displacements in the actual structure. Therefore, a number of self-equilibrating 
stress systems is introduced, represented by the matrix b,, and the magnitudes of 
these systems X are determined from the generalised compatibility equations in con- 
junction with the basic stress system and the flexibility matrix for the unassembled 
elements of the structure. The total number of the self-equilibrating stress systems 
is equal to the order of redundancy of the complete structure. The correct stress 
distribution §S, satisfying both the conditions of equilibrium and compatibility, is 
then obtained from the sum of the basic and the self-equilibrating stress 
distributions, namely 


S=b,+b,X. 


The main advantage of the method developed in the present paper lies in the 
fact that the self-equilibrating stress systems only affect the structure locally, which 
means that their matrix representation is very simple, and the effort in calculating 
these systems is not excessive. Furthermore, the mixing between systems is kept 
to a minimum and the resulting equations for the solution of the self-equilibrating 
stress systems are consequently well-conditioned. 


The analysis of shell structures with cut-outs follows here a method developed 
first in Ref. 9. Fictitious members, of arbitrary thicknesses and cross-sectional 
areas, are inserted into the cut-outs, and initial strains are impressed on these 
additional members of such magnitude that their total stresses due to both loads 
and initial strains are zero. This means that they are effectively non-existent. The 
introduction of these fictitious members ensures the retention of the uniform pattern 
of various matrices which is important for computational reasons. Furthermore, it 
should be emphasised that this approach gives the ideal method for finding the 
redistribution of internal loads due to the subsequent introduction of cut-outs in 
the structure without having to repeat all the computations. 


In considering frame flexibilities it has been assumed that the depth of the 
cross section of frames is small in comparison with the minimum radius of curva- 
ture of the shell, so that stress distributions in frames could be calculated on the 
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basis of the elementary theory of bending of beams. For simplicity, only singly- 
connected frames and shells have been considered. No difficulties, however, arise 
in cases of multiply-connected shells, except that the number of redundancies must 
be increased accordingly, e.g. double-bubble fuselage structures. The case in which 
some or all frames are treated as infinitely rigid—the common assumption of rigid 
diaphragms—is particularly simple and follows directly from the general analysis 
by equating appropriate flexibility coefficients to zero. The analysis developed 
here is quite general and is capable of dealing with stress distributions due to 
thermal loading; it may also be applied, with slight modifications, to problems 
involving diffusion of loads in curved panels with either flexible or rigid frames. 


The method of this paper has been used to detern: ie stress distributions in 
various shell structures reinforced by flexible frames and having cut-outs. In the 
Appendix a numerical example illustrating practical application of the present 
matrix method is presented; another numerical example may be found in Ref. 10. 


NOTATION 


Note: Bold face type is used to denote matrices. 


x,,y, rectangular co-ordinates defining r‘" point on the skin boundary 
x,y,’ rectangular co-ordinates defining r point on the frame neutral 
axis 
Ax. Ay matrices defined by equations (19) and (20) 
AX, Xray 
AY,=Yr-1— Yr 
As, length of r* segment 
I; iz, length of bay between frames i and (i+ 1) 
h_ distance defined in Fig. 2 


6. angle between the y-axis and an outward drawn normal to the skin 


wall at point r 
A,,A, areas defined in Fig. 3 


A,, swept area as radius vector from O moves over the skin wall from 
p tog 
(EI), flexural rigidity of the r segment 


r 


Young’s modulus 
I moment of inertia of the frame cross section 
B, cross-sectional area of r” stringer 
At cross-sectional area of r™ segment 
Aj,’ cross-sectional area of r segment effective in shear 
® skin panel area 
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panel thickness 
shear force in x- and y-directions respectively 


bending moments about y- and x-axis respectively (for sign 
convention see Fig. 2) 


torque about point F (see Fig. 2) suf 
applied loading 

frame loading (basic system) 

horizontal loads (frames) 

vertical loads (frames) 

concentrated torques (frames) 

stringer loads at a and b respectively (tension positive) 

direct load on r stringer (tension positive) 


shear flows due to applied shear load (basic system) 


shear flows due to applied torque (basic system) 2. 
resultant shear flows (basic system) 

resultant shear flow applied to a frame (basic system) ' fo 
shear flow between stringers r and (r+ 1) pCa 


distributed pressure loading 
bending moments on r" segment 
direct load on r segment 

shear load on r” segment 


distribution matrix for basic stress system 


matrix defined by equation (16) 


distribution matrix for self-equilibrating stress systems 
matrix defined by equation (45) 

distribution matrix for frame redundancies 

matrix defined by equation (43) 

self-equilibrating load systems 

frame redundancies 

relative generalised displacements 

initial strains 

Stress distribution in equilibrium with applied unit loads 


number of stringers in i" bay, continuous across the frame between 
bays i and (i+ 1) 


number of segments 


The Aeronautical Quarterly 


Ty 
R 
H, 
V, 
T, 
P., P, 
P, 
b, 
by 
bir 
Ri 
x 
| 
364 


MATRIX ANALYSIS OF SHELL STRUCTURES 
number of flexible frames 
k number of bays 


Or sign C,,C, constants defined by equations (35) and (41) 


Suffixes 
st stringers 
sk skins 
f frames 
Ss symmetric loading 
a anti-symmetric loading 
ext external load 


int internal load 


2 Structural Idealisation 


The idealised structure used here is based on a three-dimensional grid system 
formed by stringers and frames (see Fig. 1). The direct longitudinal stresses are 
carried by stringers only, while frames are capable of resisting bending moments, 
direct and shear loads, all in the plane of the frames. Bending stiffness of skin and 


(i- i) FRAME 


en 


Ficure 1. Idealised structure. 
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stringers and the warping stiffness of frames are assumed to be negligible. The 
shear flow in each skin panel is assumed to be constant and, therefore, direct loads 
in stringers between any two adjacent frames are varying linearly. Frames are 
divided into segments with at least one segment between two adjacent stringers. 
Only cylindrical shell structures of arbitrary section are considered, but it may be 
expected that the present analysis will be sufficiently accurate, without any modifi- 
cations, for shell structures with small degrees of longitudinal taper. Poisson's 
effect and plasticity are excluded, but some account of the effects due to plasticity 
may be taken, by using a method of successive approximations". 


In a continuous structure the exact evaluation of the effective areas concentra- 
ted along the grid lines is only possible when the complete distribution of direct 
Stresses is known. It is, therefore, necessary to anticipate the form of stress 
distribution in order to estimate the effective stringer areas. These effective areas 
may be taken, at least as a first approximation, to be equal to the actual stringer 
areas lying along the grid line together with half of the cross-sectional area of the 
skin material extending to the adjacent grid lines. It is always possible, should it 
prove necessary, to modify the original effective areas on the basis of a known 
(first-order approximation) stress distribution and then to re-calculate the stresses. 


3. Degree of Redundancy of the Structure 


The degree of redundancy of a singly-connected shell structure with flexible 
frames, free at one end and either fully built-in at the root section or with certain 
prescribed displacements there at all stringers, is given in Table I, where n; is the 
number of longitudinal stiffeners (stringers) which are continuous across the frame 
between bays i and (i+1), N is the number of flexible frames and k is the number 
of bays. 


If certain stringers are not held at the root section, the number of redundancies 
reduces accordingly. The numbers in the square brackets in Table I can vary from 
bay to bay, since some of the stringers can be interrupted at certain stations. If the 
root section is free, the general formulae can still be used provided that the last bay 
is not included under the summation sign, i.e. the suffix i extends from 1 to (k — 1). 


TABLE I 
Type of Loading and Structure | 


Shell Frames | Total 
| 
General case 3N 3N + S[n;—3] 
| | i=1 
Shell section with an axis of | 
symmetry | 
k k 
(a) Symmetrical loading S[(n,/2)—2] | 2N 2N + S[(n,;/2)—2] 
i=1 i=1 
k k 
(b) Anti-symmetrical loading S{(n,/2)—1] N N+ 3{(n,/2)-1 


i=1 i=1 
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iM FRAME 


£7 


= Mx 
a= h 
FicgurE 2. Basic stress system; FiGURE 3. Basic stress system; 
sign convention for applied loading. stringer loads and skin shear flows. 


4. Basic Stress System: b, Matrix 


4.1. SKINS AND STRINGERS 


It should be emphasised that the basic stress system, being a synthetic stress 
system, need not be associated with any physical cuts introduced in the structure 
to make it statically determinate. Since any stress system in equilibrium with the 
applied loading can be used for the b, matrix, it is preferable to select the simplest 
possible system. To illustrate a simple method of constructing the basic stress 
system, consider only bending moment M,, shear load S, and torque T; applied to 
the shell cross section at the (i+1)" frame. Select a pair of stringers, a and b, 
distance h apart, such that the straight line joining them is parallel to the y-axis 
(see Figs. 2 and 3). The bending moment M, is then assumed to be reacted by two 
stringer loads given by 


P,=—-—M,/h at (i+ 1)" frame (1) 
= —(M,+S,l;, i4,)/h at i frame 


while all other stringers remain unloaded, i.e. P,=0 for r#a or b. The shear 
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flows q,; and q.s arising from the shear load S, are calculated from 


qis= ( i+ i.)h (3) 
SyA, 
and Gos (A,+A,)h,’ ‘ (4) 


where A, and A, represent the swept areas as the radius vector from the origin F 
moves in the clockwise direction over the skin wall from a to b and from b to a, 
respectively. Equations (1) to (4) ensure that this simple basic load system is in 
equilibrium with S, and M,. Similar expressions can be used for shear S, and 
moment M,. If, in addition, there is torque 7; about the point F, there will be a 
constant shear flow round the section 


The expression (5) is the well-known equation for the Bredt-Batho shear flow. 
Combining equations (3) and (4) with equation (5), 


A+ 
Trl2-S,A,/h) 


The formulation of the b, matrix for the stringers and skins is extremely 
simple. It requires only a column matrix, by sen, denoting stringer loads and shear 
flows calculated from equations (1), (2), (6) and (7). For convenience the b, sue 
matrix can be written as 


b, st 
b, shell — ° ‘ (8) 


b, sk 


where the sub-matrices b, . and b, . denote stringer loads and skin shear flows 
respectively. The number of rows in the sub-matrix b, « is equal to twice the 
number of stringers in the complete structure, while the number of rows in the b, « 
sub-matrix is equal to the number of skin panels. 


Any open section is treated as a closed section with the missing panels filled 
with skin panels of similar thickness to the neighbouring panels. The effects of 
cut-outs are then obtained separately from the analysis of closed sections. The 
same method can be used for missing stringers or parts of a frame caused by the 
presence of a cut-out (see Section 8). 
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_Ficure 4. Basic stress system; frame geometry and loading. 


4.2. FRAMES 


Frames are assumed to consist of a series of segments of lengths As, measured 
along the neutral axis of the frame (see Fig. 4) and, to simplify the analysis, the 
number of segments on each frame will be taken as equal to the number of stringers 
at the section, ie. m=n. The internal loading for the basic stress system in a 
frame is described by bending moments M,,, direct loads N,,, and shear loads S,,. 
The external loading on a frame is assumed to consist of concentrated horizontal 
and vertical loads, H, and V,, and concentrated couples, 7,, all applied at the skin 
boundary at the end of each segment. Any distributed pressure or shear (shear flow) 
is replaced by corresponding concentrated loads by means of a standard matrix 
transformation. 


For the purpose of calculating the basic stress system in a frame, a “cut” is 
introduced at the first segment (see Fig. 4), the position of which is selected 
arbitrarily. For sections with an axis of symmetry, the obvious choice is the inter- 
section of this axis with the neutral axis of the frame cross section. The bending 
moment in the frame is assumed to vary linearly between the two ends of any 
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FiGurE 5. Segment of frame looking in direction zO. 


segment, while the direct and shear loads remain constant. The bending moments 


on the r™ segment are given by 


a 
M,,= 
M.,,. b 


where M,, «= — x) +T;} 


j=0 


r 


and M.,, = (y,’ — y;) V; i, Xj) + i} 


(10) 


(11) 


are the bending moments at the ends of the r segment due to all loads applied to 
the frame either externally or internally (by shear flows). The direct and shear 


loads in the frame for the basic stress system are given by 


r 


—1 r-1 
N,,-= —cos H;+sin6, = V; 


i=0 
r-l 
and S,-= —sin 6, H;—cos 6, V, 


i=0 j=0 


Now all loads applied to the frame can be written as a column matrix 
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Equations (9), (12) and (13) represent a column matrix b,; which can be calculated 
from 


where b,; is given by 
(y,’—y,) 0 00 0 0 00 0 (100 
0 00 0 00 0 
(,/-y,) ,/-y,) 0 0 0 | @-x,) 0 0 o 210 
0 0 ; 0 0 
0 0 0 0! 
0 0 | 
—cos 6, 0 0 0 0 sin 6, 0 0 0 0 
cos cos#, 0 0 0 sin 6 sin#, 0 0 0 
0 0 0 
—cos 6, —cos 4, —cos 6, sin 4, sin 4), : sin 6, 
~sin 6, 0 00 0 ~cos 6, 0 00 
—sin 6 —sin§, 0 0 0 ' —cos 6, —cos9, 0 0 0 : 0 
0 0 0 0 
—sin 6, —sin 6, —sin 6, —cos 4, —cos 6, —cos 6, 


In the actual structure, with several frames. the corresponding b,; matrix becomes 


where the suffix i refers to the i* frame and the b,;, and R;; matrices are calculated 
from equations (16) and (14) respectively. 


The loading matrix R,; includes all loads on the frame applied either externally 
or internally, the latter arising from shear flows in skin panels reacting overall 
shears and torque on the frame. If the external loading is due to pressure giving 
tise to a distributed normal loading, p,, the following transformation matrix can 
be used for calculating concentrated external loads H,ex: and V,ex: (see Fig. 5). 
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H, ext 
V. ext Ax 
pressures at 
0 0 0 0 | 
Ax, dds O 0 0 0 
0 Ax, 0 0 0 
Ax= 0 0 0 0 (19) 
0 0 0 0 
0 0 0 0 Axn 
0 0 0 0 | 
Ay, Ay, 0 0 0 0 
0 Ay, 0 0 0 
AY=! 9 oO . , 0 0 (20) 
0 0 . 0 
0 0 Aya 


The internal loads are calculated from the differences of shear flows on either side 
of the frame considered. This difference for the i frame is given by 


Gor =[b, sx] i-.,1— [Dy sk] i, (21) 


where [b, .«xJi-,,; and «xJ}s 4, denote shear flows in the and (i+1)" 
bays respectively. The shear flows q,: applied to the i" frame are then transformed 
into. concentrated loads H, int and V, in by means of equation (22), 


H, int 1 Ax 
V. int Ay 


shear flows 


The total loading applied to a frame then becomes 
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If the frame has an axis of symmetry, only one half of the structure need be 
considered. A useful check may be used at this stage for the overall equilibrium 
conditions on any frame, namely 


V.=0 28) 
=0 


{T,-—H, (y-— yr) + V, (x, - xr)} =0 


r=0 


where H,, V, and T, are obtained from equations (23). However, equations (24) 
require values of H,, V, and T,, which were not used in equations (23). 


4.3. COMPLETE b, MATRIX 


The complete b, matrix is obtained by combining equations (8) and (17). 


b, st 
by: 


It is sometimes preferable to determine the stress distribution in the shell 
structure separately for each externally applied load R,, R., ...., etc. In 
such cases the distribution matrix b, will consist of a number of column matrices, 
one corresponding to each external load; each column matrix represents the basic 
stress distribution due to a unit applied load and is determined from equation (25). 
Hence, for the complete loading, 


basic stress distribution=b, R. . ; . (25a) 


Where R= {R,, R.,. . . .} is the externally applied loading. 


5. Self-Equilibrating Stress Systems: b, Matrix 
5.1. SKINS AND STRINGERS 


(a) General Case 

In general the minimum number of stringers necessary to form a self- 
equilibrating load system is five in an unconnected (open) circuit, other than a 
Straight line, representing skins and stringers in a shell structure. The self- 
equilibrating load systems, which will be denoted as X-systems, can be obtained 
by considering five adjacent stringers at a time; the independent systems are 
obtained by advancing the first stringer in any group by one (see Fig. 6). The 
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Figure 6. Self-equilibrating stress system; general case (section between frames i and (i+1) 
looking in direction zO). 


number of linearly indpendent systems depends on the degree of redundancy of the 
stringer-skin circuit. Thus, in an m-stringer closed section (m—3) such systems can 
be constructed. The four-stringer section is a special case for which a different 


self-equilibrating system must be adapted. This special system may be obtained 
from the eigenload system of the four-boom tube theory used in Ref. 12. 


For simplicity, suffixes 1, 2, 3,.. . ., etc. will be used to denote consecutive 
stringer positions r,r+1,r+2,... ., etc., where r determines the position of the 
first stringer in any given group of self-equilibrating loads. A general self- 
equilibrating load system must have zero total end load and zero moment and 
torque resultants. Hence, in any five-stringer system the stringer loads must 
satisfy the following relationships : — 


P, +P, +P, +P, +P, =0 
+P.x, +P,x, +P.x, +P;x,=0 
+P.y. +P.y, +P,;y;=0 
=0, 


where A,, represents the swept area as the radius vector from the origin O moves 
over the skin wall from p to q. The last equation in (26) represents zero torque 
about O, expressed ir terms of the stringer loads P,. The zero shear resultants in 
the x- and y-directions are ensured automatically by the first three equations in (26). 
Equations (26) leave the magnitude of the stringer loads undetermined and, there- 
fore, some convenient relationship between the stringer loads must be introduced 
in order to assign numerical values for the b, matrix representing the self-equilibrat- 
ing X-systems. 
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For such a relationship, assume that 
. 
and hence the equations for determining stringer loads become 


“TT Pp. 


P 
P 
P 
P 
P 


0 


The calculation of the corresponding shear flows in the skin between the stringers 
presents no difficulty. The equation of equilibrium in the z-direction is applied at 
each stringer in turn, assuming constant shear flow in any given bay, so that, for 
the shear flows in the skin panels between the i” and (i+ 1)" frames, 


—(P, +P.)/[;, 
—(P, + ix, 


| 


where /; ;,, is the length of the (i+1)"* bay. The shear flows in the i bay are 
obtained by multiplying equations (29) by the factor (— 8), where B=/;, 44,/li-1,i. A 
pictorial representation of one such self-equilibrating system for the general case is 
shown in Fig. 7. 


The self-equilibrating shear flows in the skin panels acting along the periphery 
of the section are balanced by self-equilibrating loads in the frames. Thus each 
X-system will introduce additional loads into the frames. It should be noticed 
that, in general, three frames are affected by one self-equilibrating load system. 


(b) Symmetrical Structure 


For a section symmetrical about, say, the y-axis, the calculation of the b, matrix 
is greatly simplified. Any given set of applied loads can be resolved into a set 
which is symmetrical with respect to the y-axis and a set which is anti-symmetrical 
with respect to that axis. A set of symmetrical applied loads will give rise only 
to a symmetrical self-equilibrating stress system; similarly, a set of anti-symmetrical 
applied loads will produce only an anti-symmetrical self-equilibrating stress system. 
For symmetrical stress systems the stringer loads on one side of the y-axis are a 
mirror image of those on the other side. For anti-symmetrical stress systems the 
signs of stringer loads on the other side must be changed. Obviously only one 
half of the structure need be considered for the matrix analysis. 
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LINEAR VARIATION OF DIRECT 


FiGurE 7. Self-equilibrating stress system; general case. 


(c) Symmetrical Self-Equilibrating Stress Systems 


Here only three stringers on one side of the axis of symmetry need be con- 
sidered at a time (see Fig. 8). Equations for zero torque about the z-axis and zero 
bending moment about the y-axis are automatically satisfied. The relevant 
equations for the self-equilibrating stringer loads then become 


P, -P, +P, =1 
+ +P;,y,;=0 


_ 

Hence 

(y.—y,) 

P,.= 

2(ys—y,) 


and the corresponding shear flows in the (+1) bay are given by 


Qies= —P,,/l;, i+1 (32) 


i+] 
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Figure 8. Self-equilibrating stress system; symmetrical loading. 


Shear flows in the i** bay are obtained, as before, by multiplying equations (32) by 
the factor (— 8). Equations (31) and (32) are then used to construct all the 
necessary (41-2) symmetrical systems. 


If the circuit (1, 2, 3) used for determining any of the symmetrical self- 
equilibrating load systems is a straight line at right angles to the axis of symmetry, 
then only two stringers on one side are used and the expressions for the stringer 
loads and the shear flow in the (i+ 1) bay take a particularly simple form: — 


P,.=1/2; P= -1/2. Gla) 


—1 


32 


and 


(d) Anti-Symmetrical Self-Equilibrating Stress Systems 


For the anti-symmetrical systems, groups of four stringers must be considered 
(see Fig. 9). Only the condition of zero bending moment about the x-axis is 
Satisfied automatically. The equations for determining the self-equilibrating 
stringer loads are 


Pode (+ Pack, +P 
=@ 


377 


(33) 
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S. 


J. 


from which 
(X,A 54 — + 
Pos +X,A,,5)/C, | 


(34) 
P;.= X24 14 +%,4.)/C, 
+X;A,2)/C, 
where C,=2 {A,3 (%,—*.)+ An, (X,—X5)}. ; (35) 


The shear flows are then calculated from the equilibrium conditions at each stringer 
in turn. This gives, for the (i+ 1)" bay, 


qica= i+1 | / 
-(P,, + i+} (36) 
Yssa= —(Pyat Pos li, i+; 


If the circuit used for any of the anti-symmetrical systems is a straight line 
parallel to the axis of symmetry, then only three stringers need be considered. The 
stringer loads and shear flows on one side of the axis of symmetry are then calcula- 
ted from equations (31) and (32). Similarly, for a straight line circuit at right 
angles to the axis of symmetry equations (31) and (32) can also be used, provided 
that the y co-ordinates are replaced by x co-ordinates. 


The self-equilibrating stress systems, as shown in Fig. 9, are not sufficient to 
formulate all the necessary (4n—1) anti-symmetrical systems. Only (4n—-3) 
systems can be constructed by taking four-stringer circuits on one side of the axis 
of symmetry. Thus two additional systems are required. To obtain these 
additional stress systems, construct a special anti-symmetrical system, consisting of 
six stringers with the central shear panel bisected by the axis of symmetry, as 
shown in Fig. 10. The only conditions to be satisfied in this case are zero moment 
about the y-axis and zero torque about the origin O. 


F 
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Figure 10. Special anti-sym- 
metrical self-equilibrating stress 
system. 


(37) 


and P Pah t+ 


where A,, is the swept area as the radius vector from O moves over the skin wall 
from o to r. The resultant end load on one side of the axis of symmetry is in 
equilibrium with the shear flow in the top shear panel (1’, 1), while the total end load 
acting on the whole section is equal to zero. Taking the arbitrary relationship 
between the stringer loads as 


it follows that - 
Me P.. |=| . . (39) 
Ay, Mes Ao; Pos 0 


Py = (X,A 93 = X3Ay2)/C: | 
—(X,Ao3 — X34o1)/C; ‘ (40) 
where C.= — Ay, (%2 + Age (41 — + As + (41) 


The shear flows for the (i+ 1)" bay are given by 


(P,.+ Poa li, i+} | 
Gioa= (Poa + P;,)/Ii, i+1 
Gosa= (Pa)/ li, i+1 
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The stress system shown in Fig. 10 and a similar system, formed by the skins and 
stringers in the bottom part of the section, provide the two additional anti. 
symmetrical systems which, together with the (4n—3) systems given by equations 
(34) and (36), form the necessary (4n— 1) anti-symmetrical systems. 


5.2. FRAMES 


(a) Loading due to Self-Equilibrating Shear Flows 

For simplicity only one frame will be considered. The self-equilibrating shear 
flows in the skin panels are transformed into concentrated loads by means of the 
following equation 


Ax 
(b, sk) i (43) 

Ay 
where the matrix (b, ..) ; represents the resultant shear flows acting on the i" frame 
due to all X-systems.* These shear flows, applied internally to the frame, are 
shown in Fig. 7 for a general self-equilibrating stress system. The bending 
moments, direct and shear loads in the frame due to the self-equilibrating shear 

flows are calculated from . 


b,;=b 


where 


— Yo) 0 0 | 0 
(y,’ — Yo) 0 0 (x, —x,’) 0 
(y,’— Yo) (y,’ y) 0 (x, (x, ~ 
Yo) (y.’ y,) 0 (x, (x, X.’) 


sin 6, 0 
sin 6, sin 6, 


—cos 6, 0 0 0 
-cos#, 0 0 
: 0 
—cos#, —cos6, . — cos 6, 


| 
| 
| 
| 
| 
| 
| 


sin 6, 
—sin 6, 0 —cosé, 
-siné, —sin 4, —cosé, 


—siné, . —cosé, 


*The matrix (b, ,,); can be obtained from equation (21), provided that b, ,, is replaced by b, ,.- 
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FicurE 11. Frame redundancies. 


THERAME 


fe) x 


(b) Loading due to Frame Redundancies 


In general, there are three frame redundancies X,;, X.;, and X;;, on each frame 
(see Fig. 11). For a symmetrical structure these redundancies can be taken at the 


| point where the y-axis intersects the neutral axis of the frame. For symmetrical 


loading X,;=0, while for anti-symmetrical loading X,,=X,;=0. The distribution 
of the bending moments, and the direct and shear loads on the frame due to these 
redundancies will be denoted by the matrix b,,. This matrix is given by 


0 0 1 
(y,’ Yo) - 1 
(y,’—Yo) ] 
(y.’ Yo) (X, Xs ] 
—cos 6, sin 6, 0 
cos 9. sin 6, 0 (46) 
— cos 6, sin 6,, 0 
sin 6, 4, 
— sin 4, — cos 4, 0 
— sin 6, — cos 6, 0 
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5.3. COMPLETE b, MATRIX 


Combining equations (44) and (46), together with the matrices b, « and b, x 
representing stringer loads and skin shear flows respectively, and noticing that the 
frame redundancies X,;, X.;, and X,; do not induce any loads in the skins and 
stringers, the complete b, matrix becomes 


b, 0 
b,= b, sk 0 . . . . . (47) 
b,; 


6. Flexibility Matrix: f-Matrix 

The flexibility matrix f relates the deformations of the individual unassembled 
elements to their corresponding forces. It may be conveniently partitioned into 
sub-matrices representing flexibilities of stringers, skins and frames : — 


0 0 
f= 0 sx 0 |. ‘ (48) 
0 0 fy 


The evaluation of the flexibility matrices for the direct stress-carrying members 
(stringers) and skin panels has already been fully discussed in Ref. 9. The flexibility 
matrix f,, is a diagonal matrix of 2x2 sub-matrices representing flexibility 
coefficients of the unassembled stringers. A general form of these sub-matrices is 


3EB 6EB 


6EB 3EB 


where / is the stringer length, B is the stringer cross-sectional area and E is Young's 
modulus. The flexibility matrix f,. is a diagonal matrix with elements ‘/Gi, 
where “P is the area of the shear panel, ¢ is the panel thickness and G is the modulus 
of rigidity. The frame flexibility matrix f,; can be further partitioned into bending, 
extensional and shear flexibilities so that 
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The f, matrix for the bending flexibilities of the frame segments is a diagonal matrix 
nd b of 2x 2 sub-matrices whose general form is 
1 sk 


that the = _ 
ins and 
(i-sy s(1-s) 
As, { “ED, ds As, | (ED, ds 
As (50-9 4 As ds 
(47) (ED), (ED, 


' where (EJ), is the flexural rigidity of the frame cross section at the r™ segment of 
' length As, and s is the non-dimensional distance along the segment. For constant 
frame section in a segment this matrix reduces to 


mbled 
d into As, As, 
3(ED.  6(ED, 
As, As |’ 
6 (EN), 3 (ED), 
(48) As a close approximation, this expression can be applied to a variable frame section 
_ by using an average value of (EJ), for any given segment. The extensional and 
shear flexibility matrices f, and f, are purely diagonal matrices consisting of terms 
nbers 
bility As, |(EA;,) and As, 
bility 
Ces is where A;, is the frame area effective in carrying direct load, while A;,’ is the effec- 
tive area in shear. Usually A;,’ can be taken as equal to the web area of the 
frame section. 
1. Solution of the Compatibility Equations, Stress Distribution and 
Defiections 
The unknown self-equilibrating stress systems and the frame redundancies are 
ing’s — calculated from the generalised compatibility equations obtained from the unit load 
/Gt, | theorem. The derivation of the compatibility equations in matrix notation has 
ulus |) been discussed in detail by Argyris® ' and, therefore, only a brief summary of the 
ling, | relevant results will be presented. It can be shown that the required equations, 
| from which the X-systems are determined, can be expressed simply as 
‘ ; . (50) 
(49) where v is a column matrix of the relative generalised displacements.* 


*Primes are used to indicate matrix transposition, and thus b,’ denotes the transpose of b.,. 
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The displacements v are calculated from 
where S=b,R+b,X . ; (82) 


and represents the true generalised loads in the structure. Hence the generalised 
compatibility equations, obtained from equations (50), (51) and (52), become 


and the true load distribution is then obtained from 
S=b,R+b,X=b,R-b, « (4) 


The true displacements r in the directions of the applied loading R are 
determined from 


where b represents any stress distribution in equilibrium with the unit loads R. The 
b, matrix can be used in place of b, but only if the basic system is not affected by 
any cut-outs or missing elements in the structure, otherwise b=b, +b,X, calculated 
for unit loads, should be used. Substituting equations (51) and (54) into 
equation (55), 


where F=b’f (b, — b, (b,’fb,)~ 'fb,) (57) 


is the flexibility matrix for the loads R. 


For other detailed theorems for the calculation of flexibility and stiffness 
matrices in structural analysis, Ref. 11 should be consulted. 


8. Cut-Outs, Initial Strains, Thermal Loading 


Any cut-outs in stiffened shell structures can be best analysed by the technique 
described in Ref. 9. All cut-outs or missing stringers are replaced by fictitious 
shear panels or stringers having similar dimensions to the surrounding structure. 
To obtain the same stress distribution in the modified structure as in the original 
structure with cut-outs or missing stringers, initial strains are imposed on the ficti- 
tious elements of such magnitude that the stresses in these elements become zero. 
If all matrices are partitioned so that matrix elements referring to the fictitious parts 
of the structure are separated from those representing the actual structure, and if 
the column matrix of the unknown initial strains is denoted by H, then the basic 
compatability equations are given by equation (58). 
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0 
b,'fb,X | 
H 


Hence the load distribution is given by 
0 
S=(b, — b, (b,’fb,)~ ‘fb,) R— b, (b,fb,)~*b,’ + G9) 
H 


The magnitudes of the unknown strains H are obtained from equation (59) by 
equating to zero the stresses in all missing elements. This results in a set of linear 
equations for the unknown strains H. 


A similar method can also be used for frames with pin-ended joints by assum- 
ing initial slopes, at the ends of the segments at the pin joint, in the matrix H and 
then determining their magnitudes from equations (59), using the condition of zero 
moment at the joint. No difficulties arise in the application of this method to 
thermal loading and to induced stresses due to the initial lack of fit. 


§. Rigid Frames 


The analysis developed here can be considerably simplified if the frames are 
| treated as infinitely rigid. All flexibility coefficients in the frame flexibility matrix 
f, are equated to zero, which implies that the b,, b,, and f matrices will simplify to 


b, st | 


b, sk 


(60) 


(62) 


It is interesting to note that, for shell structures with rigid frames, the present 
method is analogous to the solution of structural problems by means of an 
orthogonal set of redundancies as used by Argyris and Dunne"? in their theory of 
tubes based on the fundamental assumption of an infinitely closely-spaced system 
of rigid diaphragms. The main difference between the analytical method of Argyris 
and Dunne and the present matrix method lies in the non-orthogonal character of 
the self-equilibrating stress systems (X-systems) used here and the finite spacing 
of frames. 
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10. Conclusions 


It has been demonstrated that the determination of the stress distribution in a 
cylindrical sheil structure of arbitrary cross section reinforced by flexible frames, 
under a given set of loads R, requires only three basic matrices b,, b,, and f, which 
can be compiled using very simple formulae. The b, and b, matrices depend only 
on the geometry of the structure, while the flexibility matrix f is, to a first approxi- 
mation, a function of geometry and elastic properties. Solutions of the generalised 
compatibility equations b,’fb,X+b,’fb,R=0 and the resulting load distribution 
S=b,R+b,X can be obtained on high speed electronic digital computers, either 
directly or, if the matrices are too large for the internal memory of the computer, 
from the solution of systems constructed by partitioning the basic matrices. 
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Appendix 


MATRIX ANALYSIS OF A CIRCULAR CYLINDER 


As a numerical example, consider a uniform cantilever circular cylinder with a 
uniform edge frame loaded by a concentrated radial load. This particular case has 
been selected to simplify the numerical calculations and to make a direct comparison 
with some well-known results for circular cylinders with flexible frames. In this type 
of problem, the stress distribution is influenced primarily by two non-dimensional 
structural parameters“*), The first parameter is defined as 


3 3 
(63) 


where R is the radius of the cylinder, ¢ is the skin thickness effective in carrying direct 
stresses, J is the moment of inertia of the cross section and L is the frame pitch. The 
second and somewhat less important parameter is given by 


Et/R\? 
. . . . . (64) 


where E is Young’s modulus, G is the shear modulus and f, is the skin thickness effective 
in carrying shear stresses. The remaining parameters are 


I 
C= AR? (extensional flexibility parameter) (65) 
f 


A, 


(shear flexibility parameter) . (66) 


(eccentricity parameter) . . . (67) 
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Here A, is the frame cross-sectional area, A, is the frame area effective in shear and e 
is the distance from the skin median line to the centroid of the frame cross section. 


Because of symmetry, only one-half of the structure need be considered in the 
matrix analysis. For the idealised structure, the direct stress-carrying material of the 
skin is assumed to be concentrated at 18 points uniformly spaced around the circum- 
ference and consequently the shear panels and frame segments subtend an angle of 20°. 
As a further simplification, it is assumed that the extensional and shear flexibilities of 
the frame are negligible and that the skin median line coincides with the centroid of 
the frame cross section, ie. C=D=n=0. In this example A and B are taken as 10* 
and 200 respectively, so that a direct comparison with the results of Ref. 13 could be 
made. The following elastic and geometric constants are used in the analysis : — 


R=20in. 38-2383 in.” E=10 x 10° Ib. /in2 
t=t,—0-048 in. As=6-98132 in. G= 4x 10° Ib./in.2 
L=5-47723 in. [=11-2174 in! B cin =0°335103 in.2=2=Rt/18 


(Note: The values of L and J have been calculated from equations (63) and (64) by 
assuming A =10* and B=200.) 


The numbering system for various structural elements is shown in Fig. 12: 
numbers 1-20 denote stringer elements, 21-29 denote shear panels, and 30-47 denote 
frame segments (bending terms only). There are 10 redundancies in all: 8 self- 
equilibrating load systems (X, to X,) and 2 frame redundancies (X, and X,,). 


The compatibility equation b,’fb, X + b,’fb,R=0 for R=1 can be written fully as 


0 fy, 
0 |{X}+ 


00: by; 


Since the extensional and shear flexibilities of the frame have been neglected, f,=f,=0 
and therefore, after matrix multiplication, the matrix equation (68) reduces to 
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FicureE 12. Cantilever cylinder subjected to a concentrated radial load: numbering of 


structural elements. 


In equation (69) the b,;, b,p and b,,; matrices now involve only bending moment terms. 
Also, because of symmetry of loading, the matrix b,p consists of only two columns 
corresponding to the frame redundancies X,,; and X,, (see Fig. 11). 


The basic matrices calculated for this problem are as follows:— 


0 
00684653 

0 

0 0206263 

0206263 
0800177 
0-800177 
1:91636 
wl = 1-91636 


0 
— 0:0684653 


—0-0125 
—0-0125 
— 0-0125 
—0-0125 
— 0-0125 
—0-0125 
— 0-0125 
—0-0125 
—0-0125 


— 


3-62646 
3-62646 
5-93052 
5-93052 
8-75681 
8-75681 
11-9708 
11-9708 
| 15-3910 
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1 


0 
0-371114 
0 
0-302534 
0 0 
0-128886 —0-5 0:275451 
0 0 
0:197466 -—0-5 0-257773 
0 0 
0-224549 -—0-5 0-242227 
0 0 
0242227 0-224549 
0 0 
0:257773 0-197466 
0 0 
0-275451 0-128886 
0 0 
0-302534 —0-5 
0 
0-371114 


1 
6-77558 
2-35313 — 5-52349 
3-60522 — 5-02902 
4-09968 — 4-70627 
4-42244 4.42244 
4-70627 — 4-09968 
5-02902 — 3-60522 
5-52349 — 2-35313 
6-77558 
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18 19 20 


f,, =0-272415 x 


= 199-158 x 10-° 


f, =0-0103728 x 
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o MATRIX METHOD 
— REF. 13 ; C=D=7=0 


a7 LNs 


tp 


-0.4 


FicurE 13. Shear flow distribution in a single-bay cantilever for A=104 and B=200: 


C=D=n=0. 


° MATRIX METHOD 
— REF 13; C=D=7=0 


e° 


1 
Ico 120 140 


-0.8 


Ficure 14, Bending moment distribution in the loaded frame of a single-bay cantilever cylinder 
for A= 10* and B=200; C=D=n=0. 
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MATRIX ANALYSIS OF SHELL STRUCTURES 


The row and column numbers in these sub-matrices are those of the complete b,, b, 
and f matrices. Equations (70), (71) and (72), when substituted into equation (69), 
result in ten simultaneous linear equations, the solution of which is 


X, = —0-196678 X,, = —0-570105 
X, = —0-393974 X, = — 0-275335 
X,= —0-605116 X 0-0121081 
X,= —0-745587 xX 0-0266577 


X,= —0-743075 0°442478 


The load distribution in the structure is then calculated from 


which gives the direct loads in the idealised stringer elements (and hence direct stresses 
in the skin), shear flows in the skin and bending moments in the edge frame. The 
shear flow distribution from equation (74) is plotted in Fig. 13, with the results of 
Ref. 13, while the frame bending moment is given in Fig. 14. It should be noted that, 
despite the use of a relatively crude idealisation for the continuous direct stress-carrying 
material, the shear flow and the frame bending moment agree closely with the values 
obtained in Ref. 13. All calculations were carried out om a desk machine, since the 
size of matrices and their simplicity did not warrant the use of an electronic digital 
computer. 


It is interesting to note that for a rigid frame, i.e. [=00, the frame bending flexibility 
f,=0, and equation (69) then reduces to 


D, Fat st Xx, fat Dost 
|+ =@ 
Dix fax Di cx xX, Di fx Dosx 


It can be further noted that in this case, because of symmetry of loading, X¥,= - X,, 
X,=-X,, etc. and thus a set of only four simultaneous equations need be solved. 
As would be expected, the resulting stress distribution in the shell 


agrees exactly with the engineers’ theory of bending distribution. 


It is not claimed that there is any particular advantage in using the present matrix 
method for uniform circular shell structures, since other relatively simple methods are 
available. There is no doubt, however, that for irregular structures, cut-out problems, 
thermal loading, load diffusion in shells and for many other related structural problems, 
the Matrix Force Method is a very powerful and simple tool for the calculation of 
Stresses and deflections. 
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Some Simple Results for Two-Dimensional 
Jet-Flap Aerofoils 


D. A. SPENCE 


(Royal Aircraft Establishment, Farnborough) 


SuMMaRY: It is shown from elementary considerations that the lift coefficient 
of a thin two-dimensional wing at zero incidence, with a narrow high-velocity 
jet of momentum coefficient C, issuing from its trailing edge at a (small) 
downward inclination 7, is given by 


(xC,}' 


and the loading on the chord line (0<x<1) by 


1/2 
AC,=2 


tx (1 -— x) 


(except in a certain neighbourhood of the trailing edge), for small values of C,. 
These formulae agree well with known measurements. Interpolation formulae 
for derivatives of C, at larger values of C,, based on earlier work, are also 
given: 


™ =2 (xC,)? +0-151 0-139 


=2r (14+ 0-151 +0-219 C,). 


1. Intreduction 


The pioneer experimental work on the jet flap recently carried out at the 
National Gas Turbine Establishment“ has stimulated a number of theoretical 
investigations into idealised wing-jet systems, including those by Woods’: * and 
by the present author. It appeared from these that, although the thrust on the wing 
due to the jet is given by a very simple formula, the lift and pressure distributions 
could only be obtained from the numerical solutions of certain integral equations. 
In the present note it is shown, however, that in the zero incidence case the lift 
coefficient and the pressure distribution can be obtained directly (i.e. without solving 
the flow problem fully) when the jet momentum coefficient C, is small. 


Received August 1957. 
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NOTATION 
U, undisturbed stream velocity 


V velocity in jet 

x distance from leading edge (chord length unity) 
8 thickness of jet 

p density 

p static pressure 

z jet deflection angle 

a incidence of forward portion of wing 

8B flap deflection angle 


fi(x) non-dimensional vorticity distributions representing the 
wing (i=1, 2, 3) 


K,= lim x'f; (x) 


f (x)=rf, (x) + af. (x) + Bfs (x) 
J flux of momentum in jet 
Q flux of mass in jet 
T thrust force on wing 
L lift force 
N_ force normal to forward portion of wing 
Ny force normal to flap 
Cy, Cx, Cur, C3=L/(4p,U,”) ete. 
C, pressure coefficient, =(Pjower— Pupper)/(4poU = 2f (x) 
t/c_ thickness/chord ratio of ellipse 


m modulus of transformation between ellipse and straight line 


2. “Thin Jet” Condition 


Davidson’s experiments’ showed that the horizontal thrust T on a jet-flapped 
wing is given approximately by 


where J is the flux of momentum in the jet at exit, independently of the angle of 
ejection, provided that this is small. At larger angles thrust losses occurred due 
to mixing. We consider here, as in Ref. 4, the system in which equation (1) is 
exactly true, namely that in which the jet becomes vanishingly thin while retaining 
a finite momentum. 
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This idealisation is indicated in Fig. 1. A thin two-dimensional wing of 
unit chord is at zero incidence in an infinite incompressible stream of density p, 
and undisturbed velocity U,. From the trailing edge of the wing a high-velocity 
jet issues with downward inclination >. The flow within the jet is regarded as 
isentropic, with a different stagnation pressure from the main stream. The jet 
boundaries are assumed to be streamlines of both internal and external flows, with 
no mixing across them. The momentum flow per second in the jet is then J, say, 
equal to pV78, where p is the density of the jet and V, 5 are its mean velocity and 
width respectively. 


We may recapitulate from Ref. 4 a proof that variations in J along the jet are 
of order (U,/V,)?, using J,, V., p, for momentum flux, mean velocity and mean 
static pressure in the jet at infinity downstream, and J, V, p for those at any other 
point. The jet is, for the sake of simplicity, regarded as incompressible, 
i.e. p=constant. 


In the absence of mixing, the mass flow per second is also constant and may 


be written 


7, =1+ 


Hence 
Also, since the total head in the jet is constant, 
4p (V?- V.7)=DPo — p. 
Thus, substituting for V-V,, 


J 


2 P) 


=1+ 


But, since static pressures are continuous across the jet boundaries, (p,— p) is also 
a pressure difference in the main stream, and hence is proportional to 4p,U,’; 
thus 
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Consider now the limit V,—» 00, the momentum remaining finite: from (7) it 
follows that in this case 


and from (2) and (3) that 


Q=J/V — 0, 0. (9) 


For the remainder of this paper the jet is assumed to be in this limiting con- 
dition, which seems the one most appropriate to thin-aerofoil theory. Woods" 
has doubted its value in cases for which C;=J/(4p,U,”) is itself small; but setting 
Q=0 is analytically permissible provided only that C, is non-zero, and physically 
one might expect the flow to be controlled primarily by the momentum of the jet, 
however small C, becomes. In fact it will be found in Fig. 4 that the observed 
variation of C;, with C, is correctly predicted for small C,. 


3. Thrust Equation 


Stratford™ has shown that at infinity downstream the jet must be parallel to 
the undisturbed stream (for otherwise the latter would have acquired an infinite 
vertical momentum component), and hence, by considering the balance between 
force and momentum on the fluid contained between two vertical planes, one far 
upstream and one far downstream of the wing, that 


independently of the angle of exit, as stated earlier. 


(More generally, Woods™ points out that when the fluid in the jet is drawn 
into the aerofoil from the oncoming stream and subsequently ejected with increased 
velocity at a mass flow rate Q per second, the thrust is given by 


When Q is non-zero, J varies in an unknown manner along the jet, and J,, its 
value at infinity, could not of course be specified in a practical case (in which, more- 
over, the thrust would also be modified by mixing of jet and free stream). This 
difficulty has been avoided here by taking Q=0, for which case, as just shown, 
J=constant.) 


4. Forces Acting on the Wing 


It is useful to consider, as in Fig. 2, a more general wing arrangement, 
incorporating a plain flap deflected at an angle 8 to the forward part of the wing, 
the latter being itself at incidence 2 to the undisturbed stream. The jet now 
emerges with deflection + from the trailing edge of the flap, i.e. at an angle (+ + 2+ §) 
to the stream direction. 
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FIGURE 2. 


The pressures on the wing may be integrated to a component S§ tangential to 
the surface at the leading edge, and components N, Ny normal to the fore and aft 
segments respectively. In addition the jet provides a reaction ' at the trailing edge 
acting Opposite to the exit direction. Thus, resolving horizontaily and vertically, 


T =S cos z—N sin z— Ny sin (z+ 8)+J cos (7+ 2+). . (il) 
L=N cos z+ Ny cos (z+ 8)+J sin (+ +2+)+S sin z. 


By equation (1), T may be replaced by J in equation (11). In addition, treating the 
angles 7, z, 8 as small, we may write 


ri 2 =y +2 
3p U, Or 


C +8 3B (13) 


or da” of 


(14) 


where the partial derivatives in equations (13) and (14), are independent of =, 2, f. 


Now, writing S=4p,U,°C,, N may be eliminated between equations (11) and 
(12) and, retaining only second-order terms in 7, z. 8, we obtain 


2C.=7°C, + Cs) +B (2 + +2r2 +2rB + 


a5) 


5. Expression for Leading Edge Suction Coefficient 


A finite tangential force S can be expected to arise from an infinity in the 
velocity distribution, of order x~', at the leading edge (x=0) such as always occurs 
in the potential theory of flow past a thin aerofoil at lift. The singularities in the 
velocity distribution near the hinge point (x=x,, say) and near the trailing edge 
(x=1), on the other hand, are insufficiently strong to produce finite tangential 
forces. At the trailing edge, for example, the velocity distribution is, strictly 
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speaking, of the Schwarz-Christoffel type, i.e. proportional to (1 — x)-"/«+?, but with 
the usual simplification of linearised boundary conditions it is found to be 
—(r/x)log.(1—x). From neither of these expressions does a finite suction force 
arise. Replacing 1 by x,, and z by f, the same remarks apply to the hinge point. 


To calculate Cs, with +, z, 8 supposed small, as before, the velocity distributions 
on upper and lower surfaces respectively may be written 


where f, (x) etc., which do not depend on +, 2, 8, may be thought of as vorticity 
distributions, and behave like x-+ for small x. 


The suction force at the leading edge is then given in magnitude by 


Cs= say, 
where K;= lim (x), jul, 2; 3 
are constants depending only on the jet strength. 


6. Equations Connecting Lift Derivatives 


The expression (17) for Cs; may now be used on the left hand side of equation 
(15). This yields an identity which must be true for all small values of 7, 2, £; 
hence, equating coefficients, we obtain the six relations 


OC, 


-C, (20) 
(21) 
(22) 
(23) 


(24) 
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Using these identities three at a time to eliminate pairs of the K;, we obtain 


=2C, 


together with a number of other relations not immediately required. In both these 
equations the coefficient of 2C,; on the right hand side is known from potential 
theory, for the case C; =0; in fact 


ac 


. . . (28) 


where the suffix 0 indicates potential-theory values, and 6, is the angular 
co-ordinate of the hinge point, defined by 


x,=4(1+cos6,). . ‘ (29) 
Equation (27) is the standard Kutta-Joukowski result, derived on the condition 
applying here of zero tangential force at the trailing edge, and equation (28) was 
first obtained by Glauert™ using linearised boundary conditions, i.e. treating the 


hinge as a point of stepwise discontinuity in downwash. Applying these results 
to equations (25) and (26) we obtain, for small values of C;, 


and (2) =46,°C,/t, 


the latter equation holding good only if the term C,? on the right in equation (26) 
is negligible in comparison with 2C; 0Cxr/0f, i.e. if 


7. Lift and Pressure Distribution at Zero Incidence 


From (30), the lift coefficient at zero incidence is given by 


Cy =2 
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as stated in the Summary. Similarly the pressure distribution at zero incidence may 
be obtained from (31), from which it is seen that, with z=8=0, 


1 


where AC, (x) is the difference in pressure coefficient between upper and lower 
surfaces. Hence 


AC, (x,)= 


Ix, Cxr ° (36) 


or, using equations (29) and (34) and omitting the suffix 1, 


provided only that the condition (32) is satisfied. In terms of x this becomes, 
by (29), 
1-x>C,/ 16. ‘ (38) 


Thus AC, (x) is given by equation (37) everywhere on 0 < x < 1, except in the 
immediate neighbourhood of the trailing edge, where a certain region is excluded, 
non-uniformly with respect to C,, by the inequality (38). Or, taking x ~ 1, and 
eliminating C,/(1—.x), it is seen that the limit beyond which values of AC,/7 near 
the trailing edge may not be calculated by (37) is given by 


2 


In fact, in the immediate neighbourhood of the trailing edge, as in the analysis 
of Glauert™ for a hinged flap, we should have 


AC, (x) = — * log. (1-2), 


as remarked earlier. 


8. Interpolation Formulae for Larger Values of C, 


In Ref. 4 the flow problem for the present model was solved more fully, by 
finding the functions f, (x), f, (x) in equation (16). These were obtained in Fourier 
series as solutions of integral equations relating vorticity to downwash along the 
aerofoil and jet. 9C,/@r, 0C;,/d2 were computed from the solutions for values of 
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FIGURE 3. 


C,; up to 10, and it has been found possible to approximate to their values with 
good accuracy by the following interpolation formulae, whose form is suggested 
by equation (25): — 


(Ss) =4rC, ‘ (41) 
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FIGURE 4. 


These expressions are related to each other by equation (25). They involve 
two disposable constants only, i.e. the numerical coefficients in equation (41), say, 
and these were chosen to give the correct values at C;=1 and C;=4. They are 
shown by full lines in Fig. 3, being indistinguishable on this scale from the exact 
curves. The rather more cumbersome interpolation formulae proposed in the 
earlier paper are shown by broken lines. 


9. Comparison with Experimental Results for Small C, 


The experiments described by Davidson have been reported in full by 
Dimmock™ and supply data to check equations (33) and (37). The model was an 
elliptic cylinder of thickness/chord ratio 4, with a jet exit as close as practicable to 
the trailing edge. According to potential theory, the lift coefficient on such an 
ellipse should be 9/8 of that on a thin aerofoil. Accordingly the observed lift coeffi- 
cients have been plotted in the reduced form C,/{(1+t/c) 7}, (t/c=4) in Fig. 4, 
which shows data obtained for 7=31-4°. (In Ref. 1 the numerator of the corres- 
ponding expression is C, +7/(t/c) C;, which differs from C, by an amount excluded 
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from the present analysis correct to order C;*.) Up to C;=0-2 the observed values 
follow the solid curve of 2 (xC,;)! very closely, but for larger values they follow the 
broken curve, which gives the exact value according to equation (41). 


To compare pressure distributions the relation 
(C,) ellipse — m? (1 4AC,)’ (43) 


is used, where m is the modulus of transformation between the ellipse and the thin 
aerofoil, and AC, is given by (37). In terms of the eccentric angle 6, this becomes, 
in the particular case considered, 


81 (sin +k)” 


(44) 


1- 


where k=7(C,/x)}!. This has been plotted for comparison with the observed 
pressure distribution for C;=0-30, r=31-4° in Fig. 5, and the agreement is extra- 
ordinarily close everywhere ahead of the abscissa x=0:94 which is shown by a 
broken line. Behind this line (1 — x) << ~C,/16 and (44) can no longer be used. 
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Correspondence 


Moments and Deflections of a Simply-Supported Beam 
Grillage 


Correspondence from J. CLARKSON, B.Sc. (Naval Construction Research Establishment) 
on the paper by J. P. ELLINGTON and H. McCALLION published in the November 1957 
issue Of The Aeronautical Quarterly (Vol. VIII, p. 360) 


I would like to ask the authors to clarify their interpretation of the term “exact” 
in connection with grillage analysis. In principle, the analysis of unplated grillages is 
an extremely simple problem, requiring merely an application of Euler-Bernoulli beam 
theory. In its simplest form, a system of simultaneous equations is formulated, either 
in the intersection point reactions (specifying compatible displacements at the inter- 
sections) or in displacements at the intersections (specifying equilibrium at the inter- 
sections). An example of such a solution is that by Lazarides who has set up equations 
in terms of the reactions. (Lazarides, T. O.: “The Design and Analysis of Openwork 
Prestressed Concrete Beam Grillages,” Civil Engineering, 47, 552, July 1952, pp. 471-473.) 
The simultaneous equations of Lazarides may be solved to any desired degree of 
accuracy, and the solution is exact to that degree of accuracy. 


For grillages having a fairly large number of intersections, the arithmetical work 
associated with the Lazarides solution would be prohibitively lengthy using conventional 
desk calculating machines. For such problems, Holman’s series solution is of consider- 
able interest (“ A Finite Series Solution for Grillages,” Aeronautical Quarterly, Vol. 
VIII, February 1957, pp. 49-57). Holman’s solution is an improvement of an earlier 
one by Klitchieff (“Beams on Elastic Supports and on Cross-Girders,” Aeronautical 
Quarterly, Vol. Il, November 1950, pp. 157-166). It is exact, and can take account of 
a large variety of loadings; it results in sets of simultaneous equations equal in number 
to the beams in one direction, but their formulation is rather complicated. 


Since Holman’s method was developed, electronic computing has become more 
readily available, and for most digital computers a standard programme exists for 
solution of large systems of simultaneous equations. Therefore, the simple but 
numerically tedious method of Lazarides is now almost always preferred. 
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CORRESPONDENCE 


It would appear that the method of the present paper may suffer from disadvantages 
similar to those of Holman’s solution, namely, fairly complicated analysis and arithmetic. 
The writer considers that there would be considerable iabour in obtaining a numerical 
solution for a typical practical loading such as a concentrated load or uniformly 
distributed pressure. Could the authors give some indication of how many terms of 
the solution would be required for these loadings, to obtain the bending moments 
correct to within one per cent? 


Authors’ Reply to Mr. Clarkson 


The authors would like to thank Mr. Clarkson for his communication. In reply 
to his query as to the meaning of an “exact solution”, the authors meant that a closed 
solution had been obtained for a particular system of loading, and that, in effect, the 
sets of simultaneous equations in the unknown nodal moments or displacements had 
been solved without recourse to numerical analysis or assumptions beyond the usual 
beam bending theory. 


Mr. Clarkson’s comments on the use of electronic computers are highly pertinent, 
but, as in the case of much of Mr. Clarkson’s own work on grillages, at the time of 
writing computing services were not readily available to the authors. 


However, in use the method given does not involve any complicated analysis, only 
a certain amount of straightforward arithmetic. In addition, if attention is restricted to, 
say, the central bending moments and deflection, then there is no more work involved 
in analysing a large grillage than a small one. The authors have not attempted to draw 
up design charts, but these are perfectly feasible and could reduce design work 
enormously. The original intention, however, was not to propose a design method, as 
such, but to provide a readily calculated check for existing methods. 


For grillages with a large number of nodes, some of the expressions given in the 
paper can be reduced slightly. The expression for the parameter A, equation (16), can 
be expanded as 


8 


so that for values of (mz/R) up to, say, 1, then A can be taken as (k,./24k,)(nz/R)*. 


Further, for fairly small values of A, the characteristic angles z and 8 from equation 
(18) may be approximated as 


2 cosh 2 


= J (6X) ~ 2+ (6d). 
2 cos B 


Replacing cosh z by (1 + 27/2) and cos £ by (1 — 8?/2), 
(6A) (nz/R)* (k,/4k,). 


This approximation is valid for values of A up to, say, 0-02, while 2 and 2 remain 
equal for values of A up to nearly 1. Other approximations may also be developed 
on the same lines. 


In reply to the final query as to the number of terms necessary to obtain a solution 
for, say, a uniformly distributed load. the authors cannot give a categorical reply, but 
probably three or four terms would suffice. 
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CORRESPONDENCE 


Some Identities on Structural Flexibility after Buckling 


Correspondence from M. Fine (Rotol Ltd.) on the paper by E. H. MANSFIELD published 
in the August 1958 issue of The Aeronautical Quarterly (Vol. IX, p. 300). 


It is sometimes possible to approximate to a buckled configuration which satisfies 
the boundary conditions and which contains disposable parameters. If the parameters 
are determined by suitable energy criteria the constraints do no work. But if the 
parameters are determined otherwise, by collocation for example, the constraints might 
do work. And then the strain energy stored in the structure is not necessarily equal to 
the work done by the external forces. 


The paper thus provides an argument for strain energy methods—and the problem 
of estimating stiffness by other methods. 
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